Index: example/se008731_01_40pcorr.ses
===================================================================
--- example/se008731_01_40pcorr.ses	(revision 2f8ee5f914873d15682adee6093e19798c1b0d1c)
+++ example/se008731_01_40pcorr.ses	(revision 2f8ee5f914873d15682adee6093e19798c1b0d1c)
@@ -0,0 +1,51 @@
+DataFileTitle 	"Polystyrene of Markus Strobl,  Full Sine, ++ only "						
+Sample 	"Polystyrene 2 um in 53% H2O, 47% D2O "						
+Settings 	"D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast "						
+Operator 	CPD 						
+Date 	do 10 jul 2014 16:37:30 						
+ScanType 	sine one element scan 						
+Thickness [cm] 	2.00E-01						
+Q_zmax [\A^-1] 	0.05						
+Q_ymax [\A^-1] 	0.05						
+  							
+spin echo length [A] 	 Depolarisation [A-2 cm-1] 	 error depol [A-2 cm-1] 	 error SEL [A] 	 wavelength [A] 	 error wavelength [A] 	 polarisation 	 error pol 
+391.56	0.0041929	0.0036894	19.578	2.11	0.1055	1.0037	0.0032974
+1564	-0.0046571	0.0038185	78.2	2.11	0.1055	0.99586	0.003386
+2735.6	-0.017007	0.0038132	136.78	2.11	0.1055	0.98497	0.0033444
+3907.9	-0.033462	0.0035068	195.39	2.11	0.1055	0.97064	0.0030309
+5080.2	-0.047483	0.0038208	254.01	2.11	0.1055	0.9586	0.0032613
+6251.8	-0.070375	0.00376	312.59	2.11	0.1055	0.93926	0.0031446
+7423.2	-0.092217	0.0037927	371.16	2.11	0.1055	0.92117	0.0031108
+8595.5	-0.10238	0.004006	429.77	2.11	0.1055	0.91287	0.0032562
+9767.7	-0.12672	0.0038534	488.39	2.11	0.1055	0.8933	0.0030651
+10940	-0.1374	0.004243	546.98	2.11	0.1055	0.88484	0.003343
+12112	-0.16072	0.0045837	605.58	2.11	0.1055	0.86666	0.0035372
+13284	-0.16623	0.0045613	664.2	2.11	0.1055	0.86242	0.0035027
+14456	-0.18468	0.0044918	722.79	2.11	0.1055	0.84837	0.0033931
+15628	-0.19143	0.0048967	781.38	2.11	0.1055	0.84328	0.0036768
+16800	-0.20029	0.0045421	840.02	2.11	0.1055	0.83666	0.0033837
+17971	-0.19798	0.0046642	898.56	2.11	0.1055	0.83838	0.0034819
+19143	-0.21442	0.0047052	957.17	2.11	0.1055	0.82619	0.0034614
+20316	-0.20885	0.0044931	1015.8	2.11	0.1055	0.8303	0.0033218
+21488	-0.21393	0.0049186	1074.4	2.11	0.1055	0.82655	0.00362
+22660	-0.20685	0.004423	1133	2.11	0.1055	0.83179	0.0032758
+23832	-0.20802	0.0046979	1191.6	2.11	0.1055	0.83092	0.0034758
+25003	-0.19848	0.0045953	1250.2	2.11	0.1055	0.838	0.0034289
+26175	-0.21117	0.0044567	1308.8	2.11	0.1055	0.82859	0.0032881
+27347	-0.21283	0.004137	1367.4	2.11	0.1055	0.82736	0.0030477
+28520	-0.2042	0.0044587	1426	2.11	0.1055	0.83375	0.0033101
+29692	-0.2112	0.0042852	1484.6	2.11	0.1055	0.82857	0.0031615
+30864	-0.20319	0.0043483	1543.2	2.11	0.1055	0.8345	0.003231
+32036	-0.20752	0.0044297	1601.8	2.11	0.1055	0.83129	0.0032788
+33207	-0.20654	0.0043188	1660.4	2.11	0.1055	0.83201	0.0031995
+34380	-0.20126	0.0046375	1719	2.11	0.1055	0.83593	0.0034518
+35551	-0.20924	0.0042871	1777.6	2.11	0.1055	0.83001	0.0031684
+36724	-0.21323	0.0045471	1836.2	2.11	0.1055	0.82707	0.0033487
+37895	-0.21324	0.0045354	1894.7	2.11	0.1055	0.82706	0.00334
+39067	-0.19905	0.0044141	1953.4	2.11	0.1055	0.83758	0.003292
+40239	-0.1991	0.0047441	2012	2.11	0.1055	0.83754	0.003538
+41411	-0.20359	0.0050136	2070.5	2.11	0.1055	0.8342	0.003724
+42583	-0.21032	0.0049474	2129.1	2.11	0.1055	0.82922	0.0036529
+43755	-0.20689	0.0048203	2187.8	2.11	0.1055	0.83176	0.00357
+44927	-0.21075	0.0052337	2246.4	2.11	0.1055	0.8289	0.0038628
+46099	-0.19956	0.0047827	2304.9	2.11	0.1055	0.8372	0.0035653
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision b397165ea71fa5026c955689596b9f7baf7bd132)
+++ sasmodels/sesans.py	(revision 94d13f170f2a6122e9e88300306fc88b93a8abe4)
@@ -14,201 +14,61 @@
 import numpy as np  # type: ignore
 from numpy import pi, exp  # type: ignore
-from scipy.special import jv as besselj
-#import direct_model.DataMixin as model
-        
-def make_q(q_max, Rmax):
-    r"""
-    Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$
-    to $q_max$. This is the integration range of the Hankel transform; bigger range and 
-    more points makes a better numerical integration.
-    Smaller q_min will increase reliable spin echo length range. 
-    Rmax is the "radius" of the largest expected object and can be set elsewhere.
-    q_max is determined by the acceptance angle of the SESANS instrument.
+from scipy.special import j0
+
+class SesansTransform(object):
     """
-    from sas.sascalc.data_util.nxsunit import Converter
+    Spin-Echo SANS transform calculator.  Similar to a resolution function,
+    the SesansTransform object takes I(q) for the set of *q_calc* values and
+    produces a transformed dataset
 
-    q_min = dq = 0.1 * 2*pi / Rmax
-    return np.arange(q_min,
-                     Converter(q_max[1])(q_max[0],
-                                         units="1/A"),
-                     dq)
-    
-def make_all_q(data):
+    *SElength* (A) is the set of spin-echo lengths in the measured data.
+
+    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin
+    echo encoding dimension (for ToF: Q of min(R) and max(lam)).
+
+    *Rmax* (A) is the maximum size sensitivity; larger radius requires more
+    computation time.
     """
-    Return a $q$ vector suitable for calculating the total scattering cross section for
-    calculating the effect of finite acceptance angles on Time of Flight SESANS instruments.
-    If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed.
-    If the instrument has a rectangular acceptance, 2 all_q vectors are needed.
-    If the instrument has a circular acceptance, 1 all_q vector is needed
-    
-    """
-    if not data.has_no_finite_acceptance:
-        return []
-    elif data.has_yz_acceptance(data):
-        # compute qx, qy
-        Qx, Qy = np.meshgrid(qx, qy)
-        return [Qx, Qy]
-    else:
-        # else only need q
-        # data.has_z_acceptance
-        return [q]
+    #: SElength from the data in the original data units; not used by transform
+    #: but the GUI uses it, so make sure that it is present.
+    q = None  # type: np.ndarray
 
-def transform(data, q_calc, Iq_calc, qmono, Iq_mono):
-    """
-    Decides which transform type is to be used, based on the experiment data file contents (header)
-    (2016-03-19: currently controlled from parameters script)
-    nqmono is the number of q vectors to be used for the detector integration
-    """
-    nqmono = len(qmono)
-    if nqmono == 0:
-        result = call_hankel(data, q_calc, Iq_calc)
-    elif nqmono == 1:
-        q = qmono[0]
-        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono)
-    else:
-        Qx, Qy = [qmono[0], qmono[1]]
-        Qx = np.reshape(Qx, nqx, nqy)
-        Qy = np.reshape(Qy, nqx, nqy)
-        Iq_mono = np.reshape(Iq_mono, nqx, nqy)
-        qx = Qx[0, :]
-        qy = Qy[:, 0]
-        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono)
+    #: q values to calculate when computing transform
+    q_calc = None  # type: np.ndarray
 
-    return result
+    # transform arrays
+    _H = None  # type: np.ndarray
+    _H0 = None # type: np.ndarray
 
-def call_hankel(data, q_calc, Iq_calc):
-    return hankel((data.x, data.x_unit),
-                  (data.lam, data.lam_unit),
-                  (data.sample.thickness,
-                   data.sample.thickness_unit),
-                  q_calc, Iq_calc)
-  
-def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono):
-    return hankel(data.x, data.lam * 1e-9,
-                  data.sample.thickness / 10,
-                  q_calc, Iq_calc)
-                  
-def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono):
-    return hankel(data.x, data.y, data.lam * 1e-9,
-                  data.sample.thickness / 10,
-                  q_calc, Iq_calc)
-                        
-def TotalScatter(model, parameters):  #Work in progress!!
-#    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda)
-    allq = np.linspace(0,4*pi/wavelength,1000)
-    allIq = 1
-    integral = allq*allIq
-    
+    def __init__(self, z, SElength, zaccept, Rmax):
+        # type: (np.ndarray, float, float) -> None
+        #import logging; logging.info("creating SESANS transform")
+        self.q = z
+        self._set_hankel(SElength, zaccept, Rmax)
 
+    def apply(self, Iq):
+        # tye: (np.ndarray) -> np.ndarray
+        G0 = np.dot(self._H0, Iq)
+        G = np.dot(self._H.T, Iq)
+        P = G - G0
+        return P
 
-def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still
-#==============================================================================
-#     2D Cosine Transform if "wavelength" is a vector
-#==============================================================================
-#allq is the q-space needed to create the total scattering cross-section
+    def _set_hankel(self, SElength, zaccept, Rmax):
+        # type: (np.ndarray, float, float) -> None
+        # Force float32 arrays, otherwise run into memory problems on some machines
+        SElength = np.asarray(SElength, dtype='float32')
 
-    Gprime = np.zeros_like(wavelength, 'd')
-    s = np.zeros_like(wavelength, 'd')
-    sd = np.zeros_like(wavelength, 'd')
-    Gprime = np.zeros_like(wavelength, 'd')
-    f = np.zeros_like(wavelength, 'd')
-    for i, wavelength_i in enumerate(wavelength):
-        z = magfield*wavelength_i
-        allq=np.linspace() #for calculating the Q-range of the  scattering power integral
-        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
-        alldq = (allq[1]-allq[0])*1e10
-        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
-        s[i]=1-exp(-sigma)
-        for j, Iqy_j, qy_j in enumerate(qy):
-            for k, Iqz_k, qz_k in enumerate(qz):
-                Iq = np.sqrt(Iqy_j^2+Iqz_k^2)
-                q = np.sqrt(qy_j^2 + qz_k^2)
-                Gintegral = Iq*cos(z*Qz_k)
-                Gprime[i] += Gintegral
-#                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i]
-#                s[i] += 1-exp(Totalscatter(modelname)*thickness)
-#                For now, work with standard 2-phase scatter
+        #Rmax = #value in text box somewhere in FitPage?
+        q_max = 2*pi / (SElength[1] - SElength[0])
+        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
+        q = np.arange(q_min, q_max, q_min, dtype='float32')
+        dq = q_min
 
+        H0 = np.float32(dq/(2*pi)) * q
 
-                sd[i] += Iq
-        f[i] = 1-s[i]+sd[i]
-        P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]
+        repq = np.tile(q, (SElength.size, 1)).T
+        repSE = np.tile(SElength, (q.size, 1))
+        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq
 
-
-
-
-def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname):
-#==============================================================================
-#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS
-#==============================================================================
-#acceptq is the q-space needed to create limited acceptance effect
-    SElength= wavelength*magfield
-    G = np.zeros_like(SElength, 'd')
-    threshold=2*pi*theta/wavelength
-    for i, SElength_i in enumerate(SElength):
-        allq=np.linspace() #for calculating the Q-range of the  scattering power integral
-        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
-        alldq = (allq[1]-allq[0])*1e10
-        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
-        s[i]=1-exp(-sigma)
-
-        dq = (q[1]-q[0])*1e10
-        a = (x<threshold)
-        acceptq = a*q
-        acceptIq = a*Iq
-
-        G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq)
-
-#        G[i]=np.sum(integral)
-
-    G *= dq*1e10*2*pi
-
-    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
-    
-def hankel(SElength, wavelength, thickness, q, Iq):
-    r"""
-    Compute the expected SESANS polarization for a given SANS pattern.
-
-    Uses the hankel transform followed by the exponential.  The values for *zz*
-    (or spin echo length, or delta), wavelength and sample thickness should
-    come from the dataset.  $q$ should be chosen such that the oscillations
-    in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$).
-
-    *SElength* [A] is the set of $z$ points at which to compute the
-    Hankel transform
-
-    *wavelength* [m]  is the wavelength of each individual point *zz*
-
-    *thickness* [cm] is the sample thickness.
-
-    *q* [A$^{-1}$] is the set of $q$ points at which the model has been
-    computed. These should be equally spaced.
-
-    *I* [cm$^{-1}$] is the value of the SANS model at *q*
-    """
-
-    from sas.sascalc.data_util.nxsunit import Converter
-    wavelength = Converter(wavelength[1])(wavelength[0],"A")
-    thickness = Converter(thickness[1])(thickness[0],"A")
-    Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters
-    SElength = Converter(SElength[1])(SElength[0],"A")
-
-    G = np.zeros_like(SElength, 'd')
-#==============================================================================
-#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS
-#==============================================================================
-    for i, SElength_i in enumerate(SElength):
-        integral = besselj(0, q*SElength_i)*Iq*q
-        G[i] = np.sum(integral)
-    G0 = np.sum(Iq*q)
-
-    # [m^-1] step size in q, needed for integration
-    dq = (q[1]-q[0])
-
-    # integration step, convert q into [m**-1] and 2 pi circle integration
-    G *= dq*2*pi
-    G0 = np.sum(Iq*q)*dq*2*np.pi
-
-    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0))
-
-    return P
+        self.q_calc = q
+        self._H, self._H0 = H, H0
