Index: sasmodels/kernel_template.c
===================================================================
--- sasmodels/kernel_template.c	(revision f734e7d992ea6f2199f478553447ff0dbe047559)
+++ sasmodels/kernel_template.c	(revision f734e7d992ea6f2199f478553447ff0dbe047559)
@@ -0,0 +1,236 @@
+// GENERATED CODE --- DO NOT EDIT ---
+// Code is produced by sasmodels.gen from sasmodels/models/MODEL.c
+
+#ifdef __OPENCL_VERSION__
+# define USE_OPENCL
+#endif
+
+// If opencl is not available, then we are compiling a C function
+// Note: if using a C++ compiler, then define kernel as extern "C"
+#ifndef USE_OPENCL
+#  ifdef __cplusplus
+     #include <cmath>
+     #if defined(_MSC_VER)
+     #define kernel extern "C" __declspec( dllexport )
+     #else
+     #define kernel extern "C"
+     #endif
+     using namespace std;
+     inline void SINCOS(double angle, double &svar, double &cvar)
+       { svar=sin(angle); cvar=cos(angle); }
+#  else
+     #include <math.h>
+     #if defined(_MSC_VER)
+     #define kernel __declspec( dllexport )
+     #else
+     #define kernel
+     #endif
+     #define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+#  endif
+     #if defined(_MSC_VER)
+     inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
+     #endif
+#  define global
+#  define local
+#  define constant const
+#  define powr(a,b) pow(a,b)
+#  define pown(a,b) pow(a,b)
+#else
+#  ifdef USE_SINCOS
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+#  endif
+#endif
+
+// Standard mathematical constants:
+//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4,
+//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2)
+// OpenCL defines M_constant_F for float constants, and nothing if double
+// is not enabled on the card, which is why these constants may be missing
+#ifndef M_PI
+#  define M_PI 3.141592653589793
+#endif
+#ifndef M_PI_2
+#  define M_PI_2 1.570796326794897
+#endif
+#ifndef M_PI_4
+#  define M_PI_4 0.7853981633974483
+#endif
+
+// Non-standard pi/180, used for converting between degrees and radians
+#ifndef M_PI_180
+#  define M_PI_180 0.017453292519943295
+#endif
+
+
+%(DEFINES)s
+
+%(SOURCES)s
+
+/*
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+*/
+
+#ifdef IQ_KERNEL_NAME
+kernel void IQ_KERNEL_NAME(
+    global const double *q,
+    global double *result,
+    const int Nq,
+#ifdef IQ_OPEN_LOOPS
+  #ifdef USE_OPENCL
+    global double *loops_g,
+  #endif
+    local double *loops,
+    const double cutoff,
+    IQ_DISPERSION_LENGTH_DECLARATIONS,
+#endif
+    IQ_FIXED_PARAMETER_DECLARATIONS
+    )
+{
+#ifdef USE_OPENCL
+  #ifdef IQ_OPEN_LOOPS
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (IQ_DISPERSION_LENGTH_SUM)*2, 0);
+  wait_group_events(1, &e);
+  #endif
+
+  int i = get_global_id(0);
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    const double qi = q[i];
+#ifdef IQ_OPEN_LOOPS
+    double ret=0.0, norm=0.0;
+  #ifdef VOLUME_PARAMETERS
+    double vol=0.0, norm_vol=0.0;
+  #endif
+    IQ_OPEN_LOOPS
+    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
+    //  const double radius = loops[2*(radius_i)];
+    //  const double radius_w = loops[2*(radius_i)+1];
+
+    const double weight = IQ_WEIGHT_PRODUCT;
+    if (weight > cutoff) {
+      const double I = Iq(qi, IQ_PARAMETERS);
+      if (I>=0.0) { // scattering cannot be negative
+        ret += weight*I;
+        norm += weight;
+      #ifdef VOLUME_PARAMETERS
+        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
+        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
+        norm_vol += vol_weight;
+      #endif
+      }
+    //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,I); }
+    }
+    IQ_CLOSE_LOOPS
+  #ifdef VOLUME_PARAMETERS
+    if (vol*norm_vol != 0.0) {
+      ret *= norm_vol/vol;
+    }
+  #endif
+    result[i] = scale*ret/norm+background;
+#else
+    result[i] = scale*Iq(qi, IQ_PARAMETERS) + background;
+#endif
+  }
+}
+#endif
+
+
+#ifdef IQXY_KERNEL_NAME
+kernel void IQXY_KERNEL_NAME(
+    global const double *qx,
+    global const double *qy,
+    global double *result,
+    const int Nq,
+#ifdef IQXY_OPEN_LOOPS
+  #ifdef USE_OPENCL
+    global double *loops_g,
+  #endif
+    local double *loops,
+    const double cutoff,
+    IQXY_DISPERSION_LENGTH_DECLARATIONS,
+#endif
+    IQXY_FIXED_PARAMETER_DECLARATIONS
+    )
+{
+#ifdef USE_OPENCL
+  #ifdef IQXY_OPEN_LOOPS
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (IQXY_DISPERSION_LENGTH_SUM)*2, 0);
+  wait_group_events(1, &e);
+  #endif
+
+  int i = get_global_id(0);
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    const double qxi = qx[i];
+    const double qyi = qy[i];
+#ifdef IQXY_OPEN_LOOPS
+    double ret=0.0, norm=0.0;
+    #ifdef VOLUME_PARAMETERS
+    double vol=0.0, norm_vol=0.0;
+    #endif
+    IQXY_OPEN_LOOPS
+    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
+    //  const double radius = loops[2*(radius_i)];
+    //  const double radius_w = loops[2*(radius_i)+1];
+
+    const double weight = IQXY_WEIGHT_PRODUCT;
+    if (weight > cutoff) {
+
+      const double I = Iqxy(qxi, qyi, IQXY_PARAMETERS);
+      if (I>=0.0) { // scattering cannot be negative
+        // TODO: use correct angle for spherical correction
+        // Definition of theta and phi are probably reversed relative to the
+        // equation which gave rise to this correction, leading to an
+        // attenuation of the pattern as theta moves through pi/2.  Either
+        // reverse the meanings of phi and theta in the forms, or use phi
+        // rather than theta in this correction.  Current code uses cos(theta)
+        // so that values match those of sasview.
+      #ifdef IQXY_HAS_THETA
+        const double spherical_correction
+          = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0);
+        ret += spherical_correction * weight * I;
+      #else
+        ret += weight * I;
+      #endif
+        norm += weight;
+      #ifdef VOLUME_PARAMETERS
+        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
+        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
+      #endif
+        norm_vol += vol_weight;
+      }
+      //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,I); }
+    }
+    IQXY_CLOSE_LOOPS
+  #ifdef VOLUME_PARAMETERS
+    if (vol*norm_vol != 0.0) {
+      ret *= norm_vol/vol;
+    }
+  #endif
+    result[i] = scale*ret/norm+background;
+#else
+    result[i] = scale*Iqxy(qxi, qyi, IQXY_PARAMETERS) + background;
+#endif
+  }
+}
+#endif
