Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision f1765a2d6bd586c5af9e9b67ea26d8f6a76b1e87)
+++ sasmodels/convert.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -282,6 +282,6 @@
         suppress_magnetism = False
         for key in pars.keys():
-            suppress_magnetism = suppress_magnetism or (pars[key] != 0)
             if key.startswith("M0:"):
+                suppress_magnetism = suppress_magnetism or (pars[key] != 0)
                 pars[key] = 0
         if suppress_magnetism:
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 9eb36329a5763c7ea05bd24c85b4810c08c2c2e2)
+++ sasmodels/generate.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -565,5 +565,5 @@
     else:
         # Call 1D model with sqrt(qx^2 + qy^2)
-        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))")
+        #warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))")
         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters)
         pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:]
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 9eb36329a5763c7ea05bd24c85b4810c08c2c2e2)
+++ sasmodels/model_test.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -201,4 +201,7 @@
                 #({}, 0.0, None),
                 #({}, (0.0, 0.0), None),
+                # test vector form
+                ({}, [0.1]*2, [None]*2),
+                ({}, [(0.1, 0.1)]*2, [None]*2),
                 # test that ER/VR will run if they exist
                 ({}, 'ER', None),
@@ -257,5 +260,5 @@
                 if yi is None:
                     # smoke test --- make sure it runs and produces a value
-                    self.assertTrue(np.isfinite(actual_yi),
+                    self.assertTrue(not np.isnan(actual_yi),
                                     'invalid f(%s): %s' % (xi, actual_yi))
                 elif np.isnan(yi):
Index: sasmodels/models/adsorbed_layer.py
===================================================================
--- sasmodels/models/adsorbed_layer.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/adsorbed_layer.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -51,5 +51,5 @@
 """
 
-from numpy import inf, sqrt, pi, exp
+from numpy import inf, pi, exp, errstate
 
 name = "adsorbed_layer"
@@ -87,8 +87,9 @@
     ##scale by 10^-2 for units conversion to cm^-1
     #inten =  1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm)
-    aa = (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell
+    with errstate(divide='ignore'):
+        aa = ((sld_shell - sld_solvent)/density_shell * adsorbed_amount) / q
     bb = q * second_moment
     #scale by 10^-2 for units conversion to cm^-1
-    inten = 6.0e-02 * pi * volfraction * aa * aa * exp(-bb * bb) / radius
+    inten = 6.0e-02 * pi * volfraction * aa**2 * exp(-bb**2) / radius
     return inten
 Iq.vectorized =  True  # Iq accepts an array of q values
Index: sasmodels/models/be_polyelectrolyte.py
===================================================================
--- sasmodels/models/be_polyelectrolyte.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/be_polyelectrolyte.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -116,17 +116,4 @@
 
 
-def Iqxy(qx, qy, *args):
-    """
-    :param qx:   Input q_x-value
-    :param qy:   Input q_y-value
-    :param args: Remaining arguments
-    :return:     2D-Intensity
-    """
-    intensity = Iq(sqrt(qx**2 + qy**2), *args)
-    return intensity
-
-Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
-
-
 demo = dict(scale=1, background=0.1,
             contrast_factor=10.0,
Index: sasmodels/models/broad_peak.py
===================================================================
--- sasmodels/models/broad_peak.py	(revision 65279d8edd12fdacb658e74fa3541b84dcfd56e9)
+++ sasmodels/models/broad_peak.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -41,5 +41,5 @@
 """
 
-from numpy import inf, sqrt
+from numpy import inf, errstate
 
 name = "broad_peak"
@@ -87,19 +87,10 @@
     :return:               Calculated intensity
     """
-
-    inten = (porod_scale / q ** porod_exp + lorentz_scale
-             / (1.0 + (abs(q - peak_pos) * lorentz_length) ** lorentz_exp))
+    z = abs(q - peak_pos) * lorentz_length
+    with errstate(divide='ignore'):
+        inten = (porod_scale / q ** porod_exp
+                 + lorentz_scale / (1 + z ** lorentz_exp))
     return inten
 Iq.vectorized = True  # Iq accepts an array of q values
-
-def Iqxy(qx, qy, *args):
-    """
-    :param qx:   Input q_x-value
-    :param qy:   Input q_y-value
-    :param args: Remaining arguments
-    :return:     2D-Intensity
-    """
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
 
 demo = dict(scale=1, background=0,
Index: sasmodels/models/core_shell_sphere.c
===================================================================
--- sasmodels/models/core_shell_sphere.c	(revision 7d4b2ae2af901813f3ff60e60aa0eae98f4dfd5e)
+++ sasmodels/models/core_shell_sphere.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -1,6 +1,4 @@
 double form_volume(double radius, double thickness);
 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
-double Iqxy(double qx, double qy, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
-
 
 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
@@ -16,9 +14,4 @@
 }
 
-double Iqxy(double qx, double qy, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
-    const double q = sqrt(qx*qx+qy*qy);
-    return Iq(q, radius, thickness, core_sld, shell_sld, solvent_sld);
-}
-
 double form_volume(double radius, double thickness)
 {
Index: sasmodels/models/correlation_length.py
===================================================================
--- sasmodels/models/correlation_length.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/correlation_length.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -34,5 +34,5 @@
 """
 
-from numpy import inf, sqrt
+from numpy import inf, errstate
 
 name = "correlation_length"
@@ -57,16 +57,10 @@
     1D calculation of the Correlation length model
     """
-    porod = porod_scale / pow(q, exponent_p)
-    lorentz = lorentz_scale / (1.0 + pow(q * cor_length, exponent_l))
+    with errstate(divide='ignore'):
+        porod = porod_scale / q**exponent_p
+        lorentz = lorentz_scale / (1.0 + (q * cor_length)**exponent_l)
     inten = porod + lorentz
     return inten
-
-def Iqxy(qx, qy, lorentz_scale, porod_scale, cor_length, exponent_p, exponent_l):
-    """
-    2D calculation of the Correlation length model
-    There is no orientation contribution.
-    """
-    q = sqrt(qx ** 2 + qy ** 2)
-    return Iq(q, lorentz_scale, porod_scale, cor_length, exponent_p, exponent_l)
+Iq.vectorized = True
 
 # parameters for demo
Index: sasmodels/models/dab.py
===================================================================
--- sasmodels/models/dab.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/dab.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -65,10 +65,4 @@
     """
 
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    //return -1.0;
-    return Iq(sqrt(qx*qx + qy*qy), length);
-    """
-
 # ER defaults to 1.0
 
Index: sasmodels/models/fractal.c
===================================================================
--- sasmodels/models/fractal.c	(revision 12dbc9060a2d2f1972cc8df5a80cd114bdf3876f)
+++ sasmodels/models/fractal.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -1,4 +1,2 @@
-double form_volume(double radius);
-
 double Iq(double q,
           double volfraction,
@@ -8,13 +6,4 @@
           double sld_block,
           double sld_solvent);
-
-double Iqxy(double qx, double qy,
-          double volfraction,
-          double radius,
-          double fractal_dim,
-          double cor_length,
-          double sld_block,
-          double sld_solvent);
-
 
 double Iq(double q,
@@ -64,20 +53,2 @@
 }
 
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double volfraction,
-          double radius,
-          double fractal_dim,
-          double cor_length,
-          double sld_block,
-          double sld_solvent)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-        volfraction, radius,
-        fractal_dim, cor_length,
-        sld_block, sld_solvent);
-
-}
-
Index: sasmodels/models/fractal_core_shell.c
===================================================================
--- sasmodels/models/fractal_core_shell.c	(revision 679430173a9e06492b08f4e0ab362871445963f6)
+++ sasmodels/models/fractal_core_shell.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -10,14 +10,4 @@
           double frac_dim,
           double cor_length);
-
-double Iqxy(double qx, double qy,
-            double radius,
-            double thickness,
-            double core_sld,
-            double shell_sld,
-            double solvent_sld,
-            double volfraction,
-            double frac_dim,
-            double cor_length);
 
 double form_volume(double radius, double thickness)
@@ -56,26 +46,2 @@
 }
 
-double Iqxy(double qx, double qy,
-            double radius,
-            double thickness,
-            double core_sld,
-            double shell_sld,
-            double solvent_sld,
-            double volfraction,
-            double frac_dim,
-            double cor_length) {
-
-    const double q = sqrt(qx*qx+qy*qy);
-    return Iq(q,
-              radius,
-              thickness,
-              core_sld,
-              shell_sld,
-              solvent_sld,
-              volfraction,
-              frac_dim,
-              cor_length);
-
-}
-
-
Index: sasmodels/models/fuzzy_sphere.py
===================================================================
--- sasmodels/models/fuzzy_sphere.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/fuzzy_sphere.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -99,10 +99,4 @@
     """
 
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    //return -1.0;
-    return Iq(sqrt(qx*qx + qy*qy), sld, sld_solvent, radius, fuzziness);
-    """
-
 def ER(radius):
     """
Index: sasmodels/models/gauss_lorentz_gel.py
===================================================================
--- sasmodels/models/gauss_lorentz_gel.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/gauss_lorentz_gel.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -43,5 +43,5 @@
 """
 
-from numpy import inf, sqrt, exp
+from numpy import inf, exp
 
 name = "gauss_lorentz_gel"
@@ -94,17 +94,4 @@
 
 
-def Iqxy(qx, qy, *args):
-    """
-    :param qx:   Input q_x-value
-    :param qy:   Input q_y-value
-    :param args: Remaining aruments
-    :return:     2-D intensity
-    """
-
-    return Iq(sqrt(qx**2 + qy**2), *args)
-
-Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
-
-
 demo = dict(scale=1, background=0.1,
             gauss_scale_factor=100.0,
Index: sasmodels/models/gaussian_peak.py
===================================================================
--- sasmodels/models/gaussian_peak.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/gaussian_peak.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -45,10 +45,4 @@
              ]
 
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
-form_volume = """
-    return 1.0;
-    """
-
 Iq = """
     double scaled_dq = (q - q0)/sigma;
@@ -56,12 +50,4 @@
     """
 
-
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    //return -1.0;
-    return Iq(sqrt(qx*qx + qy*qy), q0, sigma);
-    """
-
-
 # VR defaults to 1.0
 
Index: sasmodels/models/guinier.py
===================================================================
--- sasmodels/models/guinier.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/guinier.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -46,8 +46,4 @@
 """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx + qy*qy), rg);
-    """
-
 # parameters for demo
 demo = dict(scale=1.0, rg=60.0)
Index: sasmodels/models/guinier_porod.py
===================================================================
--- sasmodels/models/guinier_porod.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/guinier_porod.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -63,5 +63,5 @@
 """
 
-from numpy import inf, sqrt, power, exp
+from numpy import inf, sqrt, exp, errstate
 
 name = "guinier_porod"
@@ -93,31 +93,21 @@
     """
     n = 3.0 - s
+    ms = 0.5*(m-s) # =(n-3+m)/2
+
+    # preallocate return value
+    iq = 0.0*q
 
     # Take care of the singular points
-    if rg <= 0.0:
-        return 0.0
-    if (n-3.0+m) <= 0.0:
-        return 0.0
+    if rg <= 0.0 or ms <= 0.0:
+        return iq
 
     # Do the calculation and return the function value
-    q1 = sqrt((n-3.0+m)*n/2.0)/rg
-    if q < q1:
-        iq = (1.0/power(q, (3.0-n)))*exp((-q*q*rg*rg)/n)
-    else:
-        iq = (1.0/power(q, m))*exp(-(n-3.0+m)/2.0)*power(((n-3.0+m)*n/2.0),
-                                                         ((n-3.0+m)/2.0))/power(rg, (n-3.0+m))
+    idx = q < sqrt(n*ms)/rg
+    with errstate(divide='ignore'):
+        iq[idx] = q[idx]**-s * exp(-(q[idx]*rg)**2/n)
+        iq[~idx] = q[~idx]**-m * (exp(-ms) * (n*ms/rg**2)**ms)
     return iq
 
-Iq.vectorized = False  # Iq accepts an array of q values
-
-def Iqxy(qx, qy, *args):
-    """
-    @param qx:   Input q_x-value
-    @param qy:   Input q_y-value
-    @param args: Remaining arguments
-    """
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-
-Iqxy.vectorized = False # Iqxy accepts an array of qx, qy values
+Iq.vectorized = True # Iq accepts an array of q values
 
 demo = dict(scale=1.5, background=0.5, rg=60, s=1.0, m=3.0)
Index: sasmodels/models/line.py
===================================================================
--- sasmodels/models/line.py	(revision 416609b683895f47c35abde0f58d67d675528739)
+++ sasmodels/models/line.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -66,5 +66,6 @@
     return Iq(qx, *args)*Iq(qy, *args)
 
-Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
+Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values
+
 
 tests = [
Index: sasmodels/models/linear_pearls.c
===================================================================
--- sasmodels/models/linear_pearls.c	(revision cf85329d4eadae1d4c59b4a1c0b6f34770d49de5)
+++ sasmodels/models/linear_pearls.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -2,11 +2,4 @@
 
 double Iq(double q,
-            double radius,
-            double edge_sep,
-            double num_pearls,
-            double pearl_sld,
-            double solvent_sld);
-
-double Iqxy(double qx, double qy,
             double radius,
             double edge_sep,
@@ -92,22 +85,2 @@
 	return result;
 }
-
-double Iqxy(double qx, double qy,
-            double radius,
-            double edge_sep,
-            double num_pearls,
-            double pearl_sld,
-            double solvent_sld)
-{
-	double q;
-	q = sqrt(qx*qx+qy*qy);
-
-	double result = linear_pearls_kernel(q,
-                    radius,
-                    edge_sep,
-                    num_pearls,
-                    pearl_sld,
-                    solvent_sld);
-
-	return result;
-}
Index: sasmodels/models/lorentz.py
===================================================================
--- sasmodels/models/lorentz.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/lorentz.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -48,8 +48,4 @@
 """
 
-Iqxy = """
-    return Iq(sqrt(qx*qx + qy*qy), cor_length);
-    """
-
 # parameters for demo
 demo = dict(scale=1.0, background=0.0, cor_length=50.0)
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision 679430173a9e06492b08f4e0ab362871445963f6)
+++ sasmodels/models/mass_fractal.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -5,10 +5,4 @@
           double mass_dim,
           double cutoff_length);
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double mass_dim,
-          double cutoff_length);
-
 
 static double _mass_fractal_kernel(double q,
@@ -53,16 +47,2 @@
            cutoff_length);
 }
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double radius,
-          double mass_dim,
-          double cutoff_length)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _mass_fractal_kernel(q,
-           radius,
-           mass_dim,
-           cutoff_length);
-}
-
Index: sasmodels/models/mass_surface_fractal.c
===================================================================
--- sasmodels/models/mass_surface_fractal.c	(revision 7ed702f2d0d4ecbcb0115970cdabf0d05ac08de0)
+++ sasmodels/models/mass_surface_fractal.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -6,11 +6,4 @@
           double cluster_rg,
           double primary_rg);
-
-double Iqxy(double qx, double qy,
-          double mass_dim,
-          double surface_dim,
-          double cluster_rg,
-          double primary_rg);
-
 
 static double _mass_surface_fractal_kernel(double q,
@@ -59,18 +52,2 @@
             primary_rg);
 }
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double mass_dim,
-          double surface_dim,
-          double cluster_rg,
-          double primary_rg)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _mass_surface_fractal_kernel(q,
-           mass_dim,
-           surface_dim,
-           cluster_rg,
-           primary_rg);
-}
-
Index: sasmodels/models/mono_gauss_coil.py
===================================================================
--- sasmodels/models/mono_gauss_coil.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/mono_gauss_coil.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -45,5 +45,5 @@
 """
 
-from numpy import inf, sqrt, exp
+from numpy import inf, exp, errstate
 
 name =  "mono_gauss_coil"
@@ -63,16 +63,11 @@
 def Iq(q, i_zero, radius_gyration):
     # pylint: disable = missing-docstring
-    z = (q * radius_gyration) * (q * radius_gyration)
-    if (q == 0).any():
-       inten = i_zero
-    else:
-       inten = i_zero * 2.0 * (exp(-z) + z - 1.0 ) / (z * z)
+    z = (q * radius_gyration)**2
+
+    with errstate(invalid='ignore'):
+        inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2
+        inten[q == 0] = i_zero
     return inten
-#Iq.vectorized = True # Iq accepts an array of q values
-
-def Iqxy(qx, qy, *args):
-    # pylint: disable = missing-docstring
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-#Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
+Iq.vectorized = True # Iq accepts an array of q values
 
 demo =  dict(scale = 1.0,
Index: sasmodels/models/multilayer_vesicle.c
===================================================================
--- sasmodels/models/multilayer_vesicle.c	(revision c6ca41e077aa705fcd195f07f143a6daaf729840)
+++ sasmodels/models/multilayer_vesicle.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -59,23 +59,2 @@
 }
 
-static
-double Iqxy(double qx, double qy,
-          double volfraction,
-          double radius,
-          double thick_shell,
-          double thick_solvent,
-          double sld_solvent,
-          double sld,
-          double n_pairs)
-{
-    double q = sqrt(qx*qx + qy*qy);
-
-    return multilayer_vesicle_kernel(q,
-           volfraction,
-           radius,
-           thick_shell,
-           thick_solvent,
-           sld_solvent,
-           sld,
-           n_pairs);
-}
Index: sasmodels/models/peak_lorentz.py
===================================================================
--- sasmodels/models/peak_lorentz.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/peak_lorentz.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -29,5 +29,5 @@
 """
 
-from numpy import inf, sqrt
+from numpy import inf
 
 name = "peak_lorentz"
@@ -59,12 +59,4 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
-def Iqxy(qx, qy, *args):
-    """
-        Return I(qx, qy)
-    """
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
-
-
 demo = dict(scale=100, background=1.0,
             peak_pos=0.05, peak_hwhm=0.005)
Index: sasmodels/models/pearl_necklace.c
===================================================================
--- sasmodels/models/pearl_necklace.c	(revision cf85329d4eadae1d4c59b4a1c0b6f34770d49de5)
+++ sasmodels/models/pearl_necklace.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -6,7 +6,4 @@
 
 double Iq(double q, double radius, double edge_separation,
-	double string_thickness, double number_of_pearls, double sld, 
-	double string_sld, double solvent_sld);
-double Iqxy(double qx, double qy, double radius, double edge_separation,
 	double string_thickness, double number_of_pearls, double sld, 
 	double string_sld, double solvent_sld);
@@ -139,11 +136,2 @@
 	return value*tot_vol;
 }
-
-double Iqxy(double qx, double qy, double radius, double edge_separation,
-	double string_thickness, double number_of_pearls, double sld, 
-	double string_sld, double solvent_sld)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return(Iq(q, radius, edge_separation, string_thickness, number_of_pearls, 
-		sld, string_sld, solvent_sld));
-}
Index: sasmodels/models/poly_gauss_coil.py
===================================================================
--- sasmodels/models/poly_gauss_coil.py	(revision 65279d8edd12fdacb658e74fa3541b84dcfd56e9)
+++ sasmodels/models/poly_gauss_coil.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -50,6 +50,5 @@
 """
 
-import numpy as np
-from numpy import inf, exp, power, sqrt
+from numpy import inf, exp, power
 
 name =  "poly_gauss_coil"
@@ -83,9 +82,4 @@
 Iq.vectorized =  True  # Iq accepts an array of q values
 
-def Iqxy(qx, qy, *args):
-    # pylint: disable = missing-docstring
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
-
 demo =  dict(scale = 1.0,
             i_zero = 70.0,
Index: sasmodels/models/polymer_excl_volume.py
===================================================================
--- sasmodels/models/polymer_excl_volume.py	(revision 416609b683895f47c35abde0f58d67d675528739)
+++ sasmodels/models/polymer_excl_volume.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -93,5 +93,5 @@
 """
 
-from numpy import inf, power, sqrt
+from numpy import inf, power, errstate
 from scipy.special import gammainc, gamma
 
@@ -120,28 +120,15 @@
     :return:          Calculated intensity
     """
-    nu = 1.0/porod_exp
-    u = q*q*rg*rg*(2.0*nu+1.0) * (2.0*nu+2.0)/6.0
-    o2nu = 1.0/(2.0*nu)
+    u = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0
+    with errstate(divide='ignore', invalid='ignore'):
+        upow = power(u, -0.5*porod_exp)
+        iq = (porod_exp*upow *
+              (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp, u) -
+               upow*gamma(porod_exp)*gammainc(porod_exp, u)))
+    iq[q <= 0] = 1.0
 
-    intensity = ((1.0/(nu*power(u, o2nu))) *
-                 (gamma(o2nu)*gammainc(o2nu, u) -
-                  1.0/power(u, o2nu) * gamma(porod_exp) *
-                  gammainc(porod_exp, u))) * (q > 0) + 1.0*(q <= 0)
-
-    return intensity
+    return iq
 
 Iq.vectorized = True  # Iq accepts an array of q values
-
-
-def Iqxy(qx, qy, *args):
-    """
-    :param qx:   Input q_x-value
-    :param qy:   Input q_y-value
-    :param args: Remaining arguments
-    :return:     2D-Intensity
-    """
-    return Iq(sqrt(qx**2 + qy**2), *args)
-
-Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
 
 
Index: sasmodels/models/polymer_micelle_kernel.c
===================================================================
--- sasmodels/models/polymer_micelle_kernel.c	(revision 1e9a108a9404075e31ddee58f1d75b34a48df2fb)
+++ sasmodels/models/polymer_micelle_kernel.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -10,18 +10,4 @@
         double d_penetration,
         double n_aggreg);
-
-
-double Iqxy(double qx, double qy,
-            double ndensity,
-            double v_core,
-            double v_corona,
-            double solvent_sld,
-            double core_sld,
-            double corona_sld,
-            double radius_core,
-            double radius_gyr,
-            double d_penetration,
-            double n_aggreg);
-
 
 static double micelle_spherical_kernel(double q,
@@ -100,30 +86,2 @@
             n_aggreg);
 }
-
-double Iqxy(double qx, double qy,
-            double ndensity,
-            double v_core,
-            double v_corona,
-            double solvent_sld,
-            double core_sld,
-            double corona_sld,
-            double radius_core,
-            double radius_gyr,
-            double d_penetration,
-            double n_aggreg)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return micelle_spherical_kernel(q,
-            ndensity,
-            v_core,
-            v_corona,
-            solvent_sld,
-            core_sld,
-            corona_sld,
-            radius_core,
-            radius_gyr,
-            d_penetration,
-            n_aggreg);
-
-}
-
Index: sasmodels/models/porod.py
===================================================================
--- sasmodels/models/porod.py	(revision 82923a6edf38edd9785480897f3e484a2eb80cf8)
+++ sasmodels/models/porod.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -26,5 +26,5 @@
 """
 
-from numpy import sqrt, power, inf, errstate
+from numpy import power, inf, errstate
 
 name = "porod"
@@ -47,14 +47,4 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
-def Iqxy(qx, qy, *args):
-    """
-    @param qx:   Input q_x-value
-    @param qy:   Input q_y-value
-    @param args: Remaining arguments
-    """
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-
-Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
-
 demo = dict(scale=1.5, background=0.5)
 
Index: sasmodels/models/power_law.py
===================================================================
--- sasmodels/models/power_law.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/power_law.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -27,5 +27,5 @@
 """
 
-from numpy import inf, sqrt
+from numpy import inf, errstate
 
 name = "power_law"
@@ -45,12 +45,8 @@
 def Iq(q, power):
     # pylint: disable=missing-docstring
-    inten = (q**-power)
-    return inten
+    with errstate(divide='ignore'):
+        iq = q**-power
+    return iq
 Iq.vectorized = True  # Iq accepts an array of q values
-
-def Iqxy(qx, qy, *args):
-    # pylint: disable=missing-docstring
-    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
-Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
 
 demo = dict(scale=1.0,
Index: sasmodels/models/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision b70cb3b65351f9b9fb36c82c3c8117a33b289fb2)
+++ sasmodels/models/pringle.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -3,12 +3,4 @@
 
 double Iq(double q,
-          double radius,
-          double thickness,
-          double alpha,
-          double beta,
-          double sld_pringle,
-          double sld_solvent);
-
-double Iqxy(double qx, double qy,
           double radius,
           double thickness,
@@ -187,21 +179,2 @@
     return 1.0e-4*form*M_PI*radius*radius*thickness;
 }
-
-double Iqxy(double qx, double qy,
-            double radius,
-            double thickness,
-            double alpha,
-            double beta,
-            double sld_pringle,
-            double sld_solvent)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-            radius,
-            thickness,
-            alpha,
-            beta,
-            sld_pringle,
-            sld_solvent);
-}
-
Index: sasmodels/models/raspberry.c
===================================================================
--- sasmodels/models/raspberry.c	(revision 616df8f1d3e93ff47223b402bbf8dd7c019c0b09)
+++ sasmodels/models/raspberry.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -2,9 +2,4 @@
 
 double Iq(double q, 
-          double sld_lg, double sld_sm, double sld_solvent,
-          double volfraction_lg, double volfraction_sm, double surf_fraction,
-          double radius_lg, double radius_sm, double penetration);
-
-double Iqxy(double qx, double qy,
           double sld_lg, double sld_sm, double sld_solvent,
           double volfraction_lg, double volfraction_sm, double surf_fraction,
@@ -88,18 +83,2 @@
     return f2;
 }
-
-
-
-double Iqxy(double qx, double qy,
-          double sld_lg, double sld_sm, double sld_solvent,
-          double volfraction_lg, double volfraction_sm, double surf_fraction,
-          double radius_lg, double radius_sm, double penetration)
-          
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-          sld_lg, sld_sm, sld_solvent,
-          volfraction_lg, volfraction_sm, surf_fraction,
-          radius_lg, radius_sm, penetration);
-
-}
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/sphere.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -78,10 +78,4 @@
     """
 
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    //return -1.0;
-    return Iq(sqrt(qx*qx + qy*qy), sld, sld_solvent, radius);
-    """
-
 def ER(radius):
     """
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ sasmodels/models/stacked_disks.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -13,16 +13,4 @@
           double layer_sld,
           double solvent_sld);
-
-double Iqxy(double qx, double qy,
-            double core_thick,
-            double layer_thick,
-            double radius,
-            double n_stacking,
-            double sigma_d,
-            double core_sld,
-            double layer_sld,
-            double solvent_sld,
-            double theta,
-            double phi);
 
 static
@@ -52,5 +40,5 @@
 	const double sinarg2 = qq*(halfheight+layer_thick)*cos(zi);
 
-    const double be1 = sas_J1c(besarg1);
+	const double be1 = sas_J1c(besarg1);
 	const double be2 = sas_J1c(besarg2);
 	const double si1 = sin(sinarg1)/sinarg1;
@@ -215,30 +203,2 @@
                     solvent_sld);
 }
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-            double core_thick,
-            double layer_thick,
-            double radius,
-            double n_stacking,
-            double sigma_d,
-            double core_sld,
-            double layer_sld,
-            double solvent_sld,
-            double theta,
-            double phi)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return stacked_disks_kernel_2d(q, qx/q, qy/q,
-                    core_thick,
-                    layer_thick,
-                    radius,
-                    n_stacking,
-                    sigma_d,
-                    core_sld,
-                    layer_sld,
-                    solvent_sld,
-                    theta,
-                    phi);
-}
-
Index: sasmodels/models/star_polymer.c
===================================================================
--- sasmodels/models/star_polymer.c	(revision 55b283e83d5d634117ba0d673c0f2121cce7ce82)
+++ sasmodels/models/star_polymer.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -2,7 +2,4 @@
 
 double Iq(double q, double radius2, double arms);
-
-double Iqxy(double qx, double qy, double radius2, double arms);
-
 
 static double _mass_fractal_kernel(double q, double radius2, double arms)
@@ -28,10 +25,2 @@
     return _mass_fractal_kernel(q, radius2, arms);
 }
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy, double radius2, double arms)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _mass_fractal_kernel(q, radius2, arms);
-}
-
Index: sasmodels/models/surface_fractal.c
===================================================================
--- sasmodels/models/surface_fractal.c	(revision 679430173a9e06492b08f4e0ab362871445963f6)
+++ sasmodels/models/surface_fractal.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -5,9 +5,4 @@
           double surface_dim,
           double cutoff_length);
-
-double Iqxy(double qx, double qy,
-            double radius,
-            double surface_dim,
-            double cutoff_length);
 
 static double _surface_fractal_kernel(double q,
@@ -49,13 +44,2 @@
     return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
 }
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-    double radius,
-    double surface_dim,
-    double cutoff_length)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _surface_fractal_kernel(q, radius, surface_dim, cutoff_length);
-}
-
Index: sasmodels/models/teubner_strey.py
===================================================================
--- sasmodels/models/teubner_strey.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/teubner_strey.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -68,19 +68,8 @@
 
 
-def form_volume(radius):
-    return 1.0
-
 def Iq(q, a2, c1, c2):
     return 1. / np.polyval([c2, c1, a2], q**2)
 Iq.vectorized = True  # Iq accepts an array of q values
 
-def Iqxy(qx, qy, a2, c1, c2):
-    return Iq(sqrt(qx**2 + qy**2), a2, c1, c2)
-Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
-
-# ER defaults to 0.0
-
-# VR defaults to 1.0
-
 demo = dict(scale=1, background=0, a2=0.1, c1=-30.0, c2=5000.0)
 tests = [[{}, 0.2, 0.145927536232]]
Index: sasmodels/models/two_lorentzian.py
===================================================================
--- sasmodels/models/two_lorentzian.py	(revision 416609b683895f47c35abde0f58d67d675528739)
+++ sasmodels/models/two_lorentzian.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -34,5 +34,5 @@
 """
 
-from numpy import inf, power, sqrt
+from numpy import inf, power
 
 name = "two_lorentzian"
@@ -92,17 +92,4 @@
 
 Iq.vectorized = True  # Iq accepts an array of q values
-
-
-def Iqxy(qx, qy, *args):
-    """
-    :param qx:   Input q_x-value
-    :param qy:   Input q_y-value
-    :param args: Remaining arguments
-    :return:     2D-Intensity
-    """
-
-    return Iq(sqrt(qx**2 + qy**2), *args)
-
-Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
 
 
Index: sasmodels/models/two_power_law.py
===================================================================
--- sasmodels/models/two_power_law.py	(revision ec45c4f5ff6e17014d968da6eb80308c2c8791f0)
+++ sasmodels/models/two_power_law.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -44,8 +44,5 @@
 """
 
-from numpy import power
-from numpy import sqrt
-from numpy import inf
-from numpy import concatenate
+from numpy import inf, power, empty, errstate
 
 name = "two_power_law"
@@ -92,28 +89,13 @@
     :return:                    Calculated intensity
     """
-
-    #Two sub vectors are created to treat crossover values
-    q_lower = q[q <= crossover]
-    q_upper = q[q > crossover]
-    coefficent_2 = coefficent_1*power(crossover, -1.0*power_1)/power(crossover, -1.0*power_2)
-    intensity_lower = coefficent_1*power(q_lower, -1.0*power_1)
-    intensity_upper = coefficent_2*power(q_upper, -1.0*power_2)
-    intensity = concatenate((intensity_lower, intensity_upper), axis=0)
-
-    return intensity
+    iq = empty(q.shape, 'd')
+    idx = (q <= crossover)
+    with errstate(divide='ignore'):
+        coefficent_2 = coefficent_1 * power(crossover, power_2 - power_1)
+        iq[idx] = coefficent_1 * power(q[idx], -power_1)
+        iq[~idx] = coefficent_2 * power(q[~idx], -power_2)
+    return iq
 
 Iq.vectorized = True  # Iq accepts an array of q values
-
-def Iqxy(qx, qy, *args):
-    """
-    :param qx:   Input q_x-value
-    :param qy:   Input q_y-value
-    :param args: Remaining arguments
-    :return:     2D-Intensity
-    """
-
-    return Iq(sqrt(qx**2 + qy**2), *args)
-
-Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
 
 demo = dict(scale=1, background=0.0,
Index: sasmodels/models/unified_power_Rg.py
===================================================================
--- sasmodels/models/unified_power_Rg.py	(revision 42356c8890c04b71169941b12a36fc0e07e2034b)
+++ sasmodels/models/unified_power_Rg.py	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -62,5 +62,5 @@
 
 import numpy as np
-from numpy import inf, exp, sqrt
+from numpy import inf, exp, sqrt, errstate
 from scipy.special import erf
 
@@ -76,14 +76,17 @@
     ilevel = int(level)
     if ilevel == 0:
-        return 1./q
+        with errstate(divide='ignore'):
+            return 1./q
 
-    result = np.zeros_like(q)
-    for i in range(ilevel):
-        exp_now = exp(-(q*rg[i])**2/3.)
-        pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i]
-        exp_next = exp(-(q*rg[i+1])**2/3.) if i < ilevel-1 else 1.
-        result += G[i]*exp_now + B[i]*exp_next*pow_now
+    with errstate(divide='ignore', invalid='ignore'):
+        result = np.empty(q.shape, 'd')
+        for i in range(ilevel):
+            exp_now = exp(-(q*rg[i])**2/3.)
+            pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i]
+            exp_next = exp(-(q*rg[i+1])**2/3.) if i < ilevel-1 else 1.
+            result += G[i]*exp_now + B[i]*exp_next*pow_now
     result[q==0] = np.sum(G[:ilevel])
     return result
+
 Iq.vectorized = True
 
Index: sasmodels/models/vesicle.c
===================================================================
--- sasmodels/models/vesicle.c	(revision 062db5a05f26e70feea0a4fd1ce809f6c457dedf)
+++ sasmodels/models/vesicle.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
@@ -2,8 +2,4 @@
 
 double Iq(double q, 
-          double sld, double sld_solvent, double volfraction,
-          double radius, double thickness);
-
-double Iqxy(double qx, double qy,
           double sld, double sld_solvent, double volfraction,
           double radius, double thickness);
@@ -49,15 +45,2 @@
     return(f2);
 }
-
-
-double Iqxy(double qx, double qy,
-          double sld, double sld_solvent, double volfraction,
-          double radius, double thickness)
-          
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-        sld, sld_solvent, volfraction,
-        radius,thickness);
-
-}
