Changes in / [cdb867f:447e9aa] in sasmodels

Files:

• example/se008731_01.ses (deleted)
• example/se008731_01_40pcorr.ses (deleted)
• example/se008731_01_40points.ses (deleted)
• example/twomicronsphere.ses (deleted)
• sasmodels/direct_model.py (modified) (2 diffs)
• sasmodels/kernelpy.py (modified) (2 diffs)
• sasmodels/resolution.py (modified) (3 diffs)
• sasmodels/sesans.py (modified) (1 diff)

sasmodels/direct_model.py

@@ -31 +31 @@
 from . import resolution2d
 from .details import make_kernel_args, dispersion_mesh
-from sas.sasgui.perspectives.fitting.fitpage import FitPage
 
 try:
@@ -194 +193 @@
         # interpret data
         if hasattr(data, 'lam'):
-        #if not FitPage.no_transform.GetValue(): #if the no_transform radio button is not active DOES NOT WORK! not active before fitting
             self.data_type = 'sesans'
         elif hasattr(data, 'qx_data'):

sasmodels/kernelpy.py

@@ -15 +15 @@
 from .generate import F64
 from .kernel import KernelModel, Kernel
-import gc
 
 try:
@@ -80 +79 @@
         """
         self.q = None
-        gc.collect()
 
 class PyKernel(Kernel):

sasmodels/resolution.py

@@ -9 +9 @@
 from numpy import sqrt, log, log10, exp, pi  # type: ignore
 import numpy as np  # type: ignore
-
-from sasmodels import sesans
-from sasmodels.sesans import SesansTransform as SesansTransform
 
 __all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D",
@@ -46 +43 @@
         raise NotImplementedError("Subclass does not define the apply function")
 
+
 class Perfect1D(Resolution):
     """
@@ -57 +55 @@
     def apply(self, theory):
         return theory
+
 
 class Pinhole1D(Resolution):

sasmodels/sesans.py

@@ -14 +14 @@
 import numpy as np  # type: ignore
 from numpy import pi, exp  # type: ignore
-from scipy.special import j0
-#from mpmath import j0 as j0
+from scipy.special import jv as besselj
+#import direct_model.DataMixin as model
         
-class SesansTransform(object):
-    #: Set of spin-echo lengths in the measured data
-    SE = None  # type: np.ndarray
-    #: Maximum acceptance of scattering vector in the spin echo encoding dimension (for ToF: Q of min(R) and max(lam))
-    zaccept = None # type: float
-    #: Maximum size sensitivity; larger radius requires more computation
-    Rmax = None  # type: float
-    #: q values to calculate when computing transform
-    q = None  # type: np.ndarray
-
-    # transform arrays
-    _H = None  # type: np.ndarray
-    _H0 = None # type: np.ndarray
-
-    def set_transform(self, SE, zaccept, Rmax):
-        if self.SE is None or len(SE) != len(self.SE) or np.any(SE != self.SE) or zaccept != self.zaccept or Rmax != self.Rmax:
-            self.SE, self.zaccept, self.Rmax = SE, zaccept, Rmax
-            self._set_q()
-            self._set_hankel()
-
-    def apply(self, Iq):
-        G0 = np.dot(self._H0, Iq)
-        G = np.dot(self._H.T, Iq)
-        P = G - G0
-        return P
-
-    def _set_q(self):
-        #q_min = dq = 0.1 * 2*pi / self.Rmax
-
-        q_max = 2*pi / (self.SE[1]-self.SE[0])
-        q_min = dq = 0.1 *2*pi / (np.size(self.SE) * self.SE[-1])
-
-        #q_min = dq = q_max / 100000
-        q=np.arange(q_min, q_max, q_min)
-        self.q = q
-        self.dq = dq
-
-    def _set_hankel(self):
-        #Rmax = #value in text box somewhere in FitPage?
-        q = self.q
-        dq = self.dq
-        SElength = self.SE
-
-        H0 = dq / (2 * pi) * q
-        q=np.array(q,dtype='float32')
-        SElength=np.array(SElength,dtype='float32')
-
-        # Using numpy tile, dtype is conserved
-        repq=np.tile(q,(SElength.size,1))
-        repSE=np.tile(SElength,(q.size,1))
-        H = dq / (2 * pi) * j0(repSE*repq.T)*repq.T
-
-        # Using numpy meshgrid - meshgrid produces float64 from float32 inputs! Problem for 32-bit OS: Memerrors!
-        #H0 = dq / (2 * pi) * q
-        #repSE, repq = np.meshgrid(SElength, q)
-        #repq=np.array(repq,dtype='float32')
-        #repSE=np.array(repSE,dtype='float32')
-        #H = dq / (2 * pi) * j0(repSE*repq)*repq
-
-        self._H, self._H0 = H, H0
-
-class SESANS1D(SesansTransform):
-    def __init__(self, data, _H0, _H, q_calc):
-        # x values of the data (Sasview requires these to be named "q")
-        self.q = data.x
-        self._H0 = _H0
-        self._H = _H
-        # Pysmear does some checks on the smearer object, these checks require the "data" object...
-        self.data=data
-        # q values of the SAS model
-        self.q_calc = q_calc # These are the MODEL's q values used by the smearer (in this case: the Hankel transform)
-    def apply(self, theory):
-        return SesansTransform.apply(self,theory)
+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 sas.sascalc.data_util.nxsunit import Converter
+
+    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):
+    """
+    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]
+
+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)
+
+    return result
+
+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 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
+
+    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
+
+
+                sd[i] += Iq
+        f[i] = 1-s[i]+sd[i]
+        P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]
+
+
+
+
+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