Index: doc/developer/calculator.rst
===================================================================
--- doc/developer/calculator.rst	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
+++ doc/developer/calculator.rst	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -7,4 +7,5 @@
 - KERNEL declares a function to be available externally
 - KERNEL_NAME is the name of the function being declared
+- MAX_PD is the maximum depth of the polydispersity loop
 - NPARS is the number of parameters in the kernel
 - PARAMETER_TABLE is the declaration of the parameters to the kernel::
@@ -29,10 +30,9 @@
         double sld_solvent
 
-- CALL_IQ(q, nq, i, pars) is the declaration of a call to the kernel::
+- CALL_IQ(q, i, pars) is the declaration of a call to the kernel::
 
     Cylinder:
 
-        #define CALL_IQ(q, nq, i, var) \
-        Iq(q[i], \
+        #define CALL_IQ(q, i, var) Iq(q[i], \
         var.length, \
         var.radius, \
@@ -42,6 +42,5 @@
     Multi-shell cylinder:
 
-        #define CALL_IQ(q, nq, i, var) \
-        Iq(q[i], \
+        #define CALL_IQ(q, i, var) Iq(q[i], \
         var.num_shells, \
         var.length, \
@@ -50,4 +49,14 @@
         var.sld_solvent)
 
+    Cylinder2D:
+
+        #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent, \
+        var.theta, \
+        var.phi)
+
 - CALL_VOLUME(var) is similar, but for calling the form volume::
 
@@ -69,7 +78,4 @@
         inline bool constrained(p1, p2, p3) { return expression; }
         #define INVALID(var) constrained(var.p1, var.p2, var.p3)
-
-- IQ_FUNC could be Iq or Iqxy
-- IQ_PARS could be q[i] or q[2*i],q[2*i+1]
 
 Our design supports a limited number of polydispersity loops, wherein
@@ -200,2 +206,16 @@
 
 TODO: cutoff
+
+For accuracy we may want to introduce Kahan summation into the integration::
+
+
+    double accumulated_error = 0.0;
+    ...
+    #if USE_KAHAN_SUMMATION
+        const double y = next - accumulated_error;
+        const double t = ret + y;
+        accumulated_error = (t - ret) - y;
+        ret = t;
+    #else
+        ret += next;
+    #endif
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
+++ sasmodels/generate.py	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -236,4 +236,6 @@
 
 TEMPLATE_ROOT = dirname(__file__)
+
+MAX_PD = 4
 
 F16 = np.dtype('float16')
@@ -420,4 +422,13 @@
     return _template_cache[filename][1]
 
+_FN_TEMPLATE = """\
+double %(name)s(%(pars)s);
+double %(name)s(%(pars)s) {
+    %(body)s
+}
+
+
+"""
+
 def _gen_fn(name, pars, body):
     """
@@ -431,18 +442,36 @@
          }
     """
-    template = """\
-double %(name)s(%(pars)s);
-double %(name)s(%(pars)s) {
-    %(body)s
-}
-
-
-"""
     par_decl = ', '.join('double ' + p for p in pars) if pars else 'void'
-    return template % {'name': name, 'body': body, 'pars': par_decl}
-
-def _gen_call_pars(name, pars):
-    name += "."
-    return ",".join(name+p for p in pars)
+    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl}
+
+def _call_pars(prefix, pars):
+    """
+    Return a list of *prefix.parameter* from parameter items.
+    """
+    prefix += "."
+    return [prefix+p.name for p in pars]
+
+_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)",
+                           flags=re.MULTILINE)
+def _have_Iqxy(sources):
+    """
+    Return true if any file defines Iqxy.
+
+    Note this is not a C parser, and so can be easily confused by
+    non-standard syntax.  Also, it will incorrectly identify the following
+    as having Iqxy::
+
+        /*
+        double Iqxy(qx, qy, ...) { ... fill this in later ... }
+        */
+
+    If you want to comment out an Iqxy function, use // on the front of the
+    line instead.
+    """
+    for code in sources:
+        if _IQXY_PATTERN.search(code):
+            return True
+    else:
+        return False
 
 def make_source(model_info):
@@ -464,24 +493,25 @@
     # for computing volume even if we allow non-disperse volume parameters.
 
-    # Load template
-    source = [load_template('kernel_header.c')]
-
-    # Load additional sources
-    source += [open(f).read() for f in model_sources(model_info)]
-
-    # Prepare defines
-    defines = []
-
-    iq_parameters = [p.name
-                     for p in model_info['parameters'][2:]  # skip scale, background
-                     if p.name in model_info['par_set']['1d']]
-    iqxy_parameters = [p.name
-                       for p in model_info['parameters'][2:]  # skip scale, background
-                       if p.name in model_info['par_set']['2d']]
-    volume_parameters = model_info['par_type']['volume']
+    # kernel_iq assumes scale and background are the first parameters;
+    # they should be first for 1d and 2d parameter lists as well.
+    assert model_info['parameters'][0].name == 'scale'
+    assert model_info['parameters'][1].name == 'background'
+
+    # Identify parameter types
+    iq_parameters = model_info['par_type']['1d'][2:]
+    iqxy_parameters = model_info['par_type']['2d'][2:]
+    vol_parameters = model_info['par_type']['volume']
+
+    # Load templates and user code
+    kernel_header = load_template('kernel_header.c')
+    kernel_code = load_template('kernel_iq.c')
+    user_code = [open(f).read() for f in model_sources(model_info)]
+
+    # Build initial sources
+    source = [kernel_header] + user_code
 
     # Generate form_volume function, etc. from body only
     if model_info['form_volume'] is not None:
-        pnames = [p.name for p in volume_parameters]
+        pnames = [p.name for p in vol_parameters]
         source.append(_gen_fn('form_volume', pnames, model_info['form_volume']))
     if model_info['Iq'] is not None:
@@ -492,51 +522,54 @@
         source.append(_gen_fn('Iqxy', pnames, model_info['Iqxy']))
 
-    # Fill in definitions for volume parameters
-    if volume_parameters:
-        deref_vol = ",".join("v."+p.name for p in volume_parameters)
-        defines.append(('CALL_VOLUME(v)', 'form_volume(%s)\n'%deref_vol))
+    # Define the parameter table
+    source.append("#define PARAMETER_TABLE \\")
+    source.append("\\\n    ".join("double %s;"%p.name
+                                   for p in model_info['parameters'][2:]))
+
+    # Define the function calls
+    if vol_parameters:
+        refs = ",".join(_call_pars("v", vol_parameters))
+        call_volume = "#define CALL_VOLUME(v) form_volume(%s)"%refs
     else:
         # Model doesn't have volume.  We could make the kernel run a little
         # faster by not using/transferring the volume normalizations, but
         # the ifdef's reduce readability more than is worthwhile.
-        defines.append(('CALL_VOLUME(v)', '0.0'))
-
-    # Fill in definitions for Iq parameters
-    defines.append(('KERNEL_NAME', model_info['name']))
-    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters)))
-    if fixed_1d:
-        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS',
-                        ', \\\n    '.join('const double %s' % p for p in fixed_1d)))
-    # Fill in definitions for Iqxy parameters
-    defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy'))
-    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters)))
-    if fixed_2d:
-        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS',
-                        ', \\\n    '.join('const double %s' % p for p in fixed_2d)))
-    if pd_2d:
-        defines.append(('IQXY_WEIGHT_PRODUCT',
-                        '*'.join(p + '_w' for p in pd_2d)))
-        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS',
-                        ', \\\n    '.join('const int N%s' % p for p in pd_2d)))
-        defines.append(('IQXY_DISPERSION_LENGTH_SUM',
-                        '+'.join('N' + p for p in pd_2d)))
-        open_loops, close_loops = build_polydispersity_loops(pd_2d)
-        defines.append(('IQXY_OPEN_LOOPS',
-                        open_loops.replace('\n', ' \\\n')))
-        defines.append(('IQXY_CLOSE_LOOPS',
-                        close_loops.replace('\n', ' \\\n')))
-    # Need to know if we have a theta parameter for Iqxy; it is not there
-    # for the magnetic sphere model, for example, which has a magnetic
-    # orientation but no shape orientation.
-    if 'theta' in pd_2d:
-        defines.append(('IQXY_HAS_THETA', '1'))
-
-    #for d in defines: print(d)
-    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
-    sources = '\n\n'.join(source)
-    return C_KERNEL_TEMPLATE % {
-        'DEFINES': defines,
-        'SOURCES': sources,
-        }
+        call_volume = "#define CALL_VOLUME(v) 0.0"
+    source.append(call_volume)
+
+    refs = ["q[i]"] + _call_pars("v", iq_parameters)
+    call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs))
+    if _have_Iqxy(user_code):
+        # Call 2D model
+        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v", iqxy_parameters)
+        call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs))
+    else:
+        # Call 1D model with 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:]
+        call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt))
+
+    # Fill in definitions for numbers of parameters
+    source.append("#define MAX_PD %s"%MAX_PD)
+    source.append("#define NPARS %d"%(len(model_info['parameters'])-2))
+
+    # TODO: allow mixed python/opencl kernels?
+
+    # define the Iq kernel
+    source.append("#define KERNEL_NAME %s_Iq"%model_info['name'])
+    source.append(call_iq)
+    source.append(kernel_code)
+    source.append("#undef CALL_IQ")
+    source.append("#undef KERNEL_NAME")
+
+    # define the Iqxy kernel from the same source with different #defines
+    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name'])
+    source.append(call_iqxy)
+    source.append(kernel_code)
+    source.append("#undef CALL_IQ")
+    source.append("#undef KERNEL_NAME")
+
+    return '\n'.join(source)
 
 def categorize_parameters(pars):
@@ -588,5 +621,5 @@
     }
     for p in pars:
-        par_type[p.type if p.type else 'other'].append(p.name)
+        par_type[p.type if p.type else 'other'].append(p)
     return  par_type
 
@@ -605,13 +638,28 @@
     pars = [Parameter(*p) for p in model_info['parameters']]
     # Fill in the derived attributes
+    par_type = collect_types(pars)
+    par_type.update(categorize_parameters(pars))
     model_info['parameters'] = pars
-    partype = categorize_parameters(pars)
     model_info['limits'] = dict((p.name, p.limits) for p in pars)
-    model_info['par_type'] = collect_types(pars)
-    model_info['par_set'] = categorize_parameters(pars)
+    model_info['par_type'] = par_type
     model_info['defaults'] = dict((p.name, p.default) for p in pars)
     if model_info.get('demo', None) is None:
         model_info['demo'] = model_info['defaults']
-    model_info['has_2d'] = partype['orientation'] or partype['magnetic']
+    model_info['has_2d'] = par_type['orientation'] or par_type['magnetic']
+
+def create_default_functions(model_info):
+    """
+    Autogenerate missing functions, such as Iqxy from Iq.
+
+    This only works for Iqxy when Iq is written in python. :func:`make_source`
+    performs a similar role for Iq written in C.
+    """
+    if model_info['Iq'] is not None and model_info['Iqxy'] is None:
+        if model_info['par_type']['1d'] != model_info['par_type']['2d']:
+            raise ValueError("Iqxy model is missing")
+        Iq = model_info['Iq']
+        def Iqxy(qx, qy, **kw):
+            return Iq(np.sqrt(qx**2 + qy**2), **kw)
+        model_info['Iqxy'] = Iqxy
 
 def make_model_info(kernel_module):
@@ -693,4 +741,5 @@
     functions = "ER VR form_volume Iq Iqxy shape sesans".split()
     model_info.update((k, getattr(kernel_module, k, None)) for k in functions)
+    create_default_functions(model_info)
     return model_info
 
Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision 2e44ac7f2ee50d35367a9432419c61b4d5ed4b56)
+++ sasmodels/kernel_header.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -1,7 +1,7 @@
 #ifdef __OPENCL_VERSION__
 # define USE_OPENCL
+#elif defined(_OPENMP)
+# define USE_OPENMP
 #endif
-
-#define USE_KAHAN_SUMMATION 0
 
 // If opencl is not available, then we are compiling a C function
@@ -16,8 +16,8 @@
          #include <float.h>
          #define kernel extern "C" __declspec( dllexport )
-         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
-         inline double fmin(double x, double y) { return x>y ? y : x; }
-         inline double fmax(double x, double y) { return x<y ? y : x; }
-         inline double isnan(double x) { return _isnan(x); }
+         static double trunc(double x) { return x>=0?floor(x):-floor(-x); }
+         static double fmin(double x, double y) { return x>y ? y : x; }
+         static double fmax(double x, double y) { return x<y ? y : x; }
+         static double isnan(double x) { return _isnan(x); }
          #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
          static double cephes_expm1(double x) {
@@ -48,5 +48,5 @@
          #define kernel extern "C"
      #endif
-     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
+     static void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
 #  else
      #include <stdio.h>
@@ -99,8 +99,6 @@
 #  define M_4PI_3 4.18879020478639
 #endif
-//inline double square(double x) { return pow(x,2.0); }
-//inline double square(double x) { return pown(x,2); }
-inline double square(double x) { return x*x; }
-inline double cube(double x) { return x*x*x; }
-inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+static double square(double x) { return x*x; }
+static double cube(double x) { return x*x*x; }
+static double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
+++ sasmodels/kernel_iq.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -12,7 +12,8 @@
 */
 
+#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
+#define _PAR_BLOCK_
 
 #define MAX_PD 4  // MAX_PD is the max number of polydisperse parameters
-#define PD_2N 16  // PD_2N is the size of the coordination step table
 
 typedef struct {
@@ -31,10 +32,11 @@
 
 typedef struct {
-    PARAMETER_DECL;
+    PARAMETER_TABLE;
 } ParameterBlock;
-
-#define FULL_KERNEL_NAME KERNEL_NAME ## _ ## IQ_FUNC
-KERNEL
-void FULL_KERNEL_NAME(
+#endif
+
+
+kernel
+void KERNEL_NAME(
     int nq,                 // number of q values
     global const ProblemDetails *problem,
@@ -42,8 +44,8 @@
     global const double *pars,
     global const double *q, // nq q values, with padding to boundary
-    global double *result,  // nq return values, again with padding
+    global double *result,  // nq+3 return values, again with padding
     const double cutoff,    // cutoff in the polydispersity weight product
     const int pd_start,     // where we are in the polydispersity loop
-    const int pd_stop,      // where we are stopping in the polydispersity loop
+    const int pd_stop       // where we are stopping in the polydispersity loop
     )
 {
@@ -52,5 +54,5 @@
   // walk the polydispersity cube.
   local ParameterBlock local_pars;  // current parameter values
-  const double *parvec = &local_pars;  // Alias named parameters with a vector
+  double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector
 
   local int offset[NPARS-2];
@@ -60,5 +62,5 @@
     // Shouldn't need to copy!!
     for (int k=0; k < NPARS; k++) {
-      parvec[k] = pars[k+2];  // skip scale and background
+      pvec[k] = pars[k+2];  // skip scale and background
     }
 
@@ -68,5 +70,5 @@
     for (int i=0; i < nq; i++) {
     {
-      const double scattering = IQ_FUNC(IQ_PARS, IQ_PARAMETERS);
+      const double scattering = CALL_IQ(q, i, local_pars);
       result[i] += pars[0]*scattering + pars[1];
     }
@@ -99,13 +101,14 @@
     norm = result[nq];
     vol = result[nq+1];
-    norm_vol = results[nq+2];
+    norm_vol = result[nq+2];
   }
 
   // Location in the polydispersity hypercube, one index per dimension.
-  local int pd_index[PD_MAX];
+  local int pd_index[MAX_PD];
 
   // Trigger the reset behaviour that happens at the end the fast loop
   // by setting the initial index >= weight vector length.
-  pd_index[0] = pd_length[0];
+  pd_index[0] = problem->pd_length[0];
+
 
   // need product of weights at every Iq calc, so keep product of
@@ -120,35 +123,37 @@
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
     // check if indices need to be updated
-    if (pd_index[0] >= pd_length[0]) {
+    if (pd_index[0] >= problem->pd_length[0]) {
 
       // RESET INDICES
-      pd_index[0] = loop_index%pd_length[0];
+      pd_index[0] = loop_index%problem->pd_length[0];
       partial_weight = 1.0;
       partial_volweight = 1.0;
       for (int k=1; k < MAX_PD; k++) {
-        pd_index[k] = (loop_index%pd_length[k])/pd_stride[k];
-        const double wi = weights[pd_offset[0]+pd_index[0]];
+        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
+        const double wi = weights[problem->pd_offset[0]+pd_index[0]];
         partial_weight *= wi;
-        if (pd_isvol[k]) partial_volweight *= wi;
+        if (problem->pd_isvol[k]) partial_volweight *= wi;
       }
       for (int k=0; k < NPARS; k++) {
-        int coord = par_coord[k];
-        int this_offset = par_offset[k];
+        int coord = problem->par_coord[k];
+        int this_offset = problem->par_offset[k];
         int block_size = 1;
         for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
           if (coord&1) {
               this_offset += block_size * pd_index[bit];
-              block_size *= pd_length[bit];
+              block_size *= problem->pd_length[bit];
           }
           coord /= 2;
         }
         offset[k] = this_offset;
-        parvec[k] = pars[this_offset];
-      }
-      weight = partial_weight * weights[pd_offset[0]+pd_index[0]]
-      if (theta_var >= 0) {
-        spherical_correction = fabs(cos(M_PI_180*parvec[theta_var]));
-      }
-      if (!fast_theta) weight *= spherical_correction;
+        pvec[k] = pars[this_offset];
+      }
+      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
+      if (problem->theta_var >= 0) {
+        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_var]));
+      }
+      if (!problem->fast_theta) {
+        weight *= spherical_correction;
+      }
 
     } else {
@@ -156,13 +161,14 @@
       // INCREMENT INDICES
       pd_index[0] += 1;
-      const double wi = weights[pd_offset[0]+pd_index[0]];
+      const double wi = weights[problem->pd_offset[0]+pd_index[0]];
       weight = partial_weight*wi;
-      if (pd_isvol[0]) vol_weight *= wi;
+      if (problem->pd_isvol[0]) vol_weight *= wi;
       for (int k=0; k < problem->fast_coord_count; k++) {
-        parvec[ fast_coord_index[k]]
-            = pars[offset[fast_coord_index[k]] + pd_index[0]];
-      }
-      if (fast_theta) weight *= fabs(cos(M_PI_180*parvec[theta_var]));
-
+        pvec[problem->fast_coord_index[k]]
+            = pars[offset[problem->fast_coord_index[k]]++];
+      }
+      if (problem->fast_theta) {
+        weight *= fabs(cos(M_PI_180*pvec[problem->theta_var]));
+      }
     }
 
@@ -180,8 +186,8 @@
       #endif
       for (int i=0; i < nq; i++) {
-      {
-        const double scattering = CALL_IQ(q, nq, i, local_pars);
+        const double scattering = CALL_IQ(q, i, local_pars);
         result[i] += weight*scattering;
       }
+    }
   }
   //Makes a normalization avialable for the next round
@@ -191,5 +197,5 @@
 
   //End of the PD loop we can normalize
-  if (pd_stop == pd_stride[MAX_PD-1]) {}
+  if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
     #ifdef USE_OPENMP
     #pragma omp parallel for
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 17bbadde9de73711781f748816b06e6956dc8cfa)
+++ sasmodels/kerneldll.py	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -61,5 +61,5 @@
 if sys.platform == 'darwin':
     #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
-    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
+    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm -fno-unknown-pragmas"
 elif os.name == 'nt':
     # call vcvarsall.bat before compiling to set path, headers, libs, etc.
@@ -80,5 +80,7 @@
         COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
         if "SAS_OPENMP" in os.environ:
-            COMPILE = COMPILE + " -fopenmp"
+            COMPILE += " -fopenmp"
+        else:
+            COMPILE += " -fWno-unknown-pragmas"
 else:
     COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 50e1e40cf14d4ded28bf96e61a3cfaf35003b899)
+++ sasmodels/models/cylinder.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -3,4 +3,6 @@
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
     double radius, double length, double theta, double phi);
+
+#define INVALID(v) (v.radius<0 || v.length<0)
 
 double form_volume(double radius, double length)
@@ -15,5 +17,4 @@
     double length)
 {
-    // TODO: return NaN if radius<0 or length<0?
     // precompute qr and qh to save time in the loop
     const double qr = q*radius;
@@ -47,5 +48,4 @@
     double phi)
 {
-    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output
 
Index: sasmodels/models/gel_fit.c
===================================================================
--- sasmodels/models/gel_fit.c	(revision 30b4ddfb88f4c0a299757def2ee628c2a3d50a42)
+++ sasmodels/models/gel_fit.c	(revision 03cac087353e872f40154bf57ef8161eb5367fd4)
@@ -1,19 +1,3 @@
-double form_volume(void);
-
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length);
-
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length);
-
-static double _gel_fit_kernel(double q,
+static double Iq(double q,
           double guinier_scale,
           double lorentzian_scale,
@@ -24,5 +8,5 @@
     // Lorentzian Term
     ////////////////////////double a(x[i]*x[i]*zeta*zeta);
-    double lorentzian_term = q*q*cor_length*cor_length;
+    double lorentzian_term = square(q*cor_length);
     lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term;
     lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) );
@@ -30,5 +14,5 @@
     // Exponential Term
     ////////////////////////double d(x[i]*x[i]*rg*rg);
-    double exp_term = q*q*gyration_radius*gyration_radius;
+    double exp_term = square(q*gyration_radius);
     exp_term = exp(-1.0*(exp_term/3.0));
 
@@ -37,39 +21,2 @@
     return result;
 }
-double form_volume(void){
-    // Unused, so free to return garbage.
-    return NAN;
-}
-
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
-{
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
-}
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
-}
-
