Changes in / [2622b3f:a792234] in sasmodels

Files:

• example/sesansfit.py (modified) (4 diffs)
• sasmodels/core.py (modified) (2 diffs)
• sasmodels/direct_model.py (modified) (6 diffs)
• sasmodels/sesans.py (modified) (3 diffs)

example/sesansfit.py

@@ -1 +1 @@
 from bumps.names import *
-from sasmodels import core, bumps_model, sesans
+from sas import core, bumps_model, sesans
 
 HAS_CONVERTER = True
@@ -8 +8 @@
     HAS_CONVERTER = False
 
+
 def get_bumps_model(model_name):
     kernel = core.load_model(model_name)
@@ -13 +14 @@
     return model
 
-def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[]):
+def sesans_fit(file, model_name, initial_vals={}, custom_params={}, param_range=[], acceptance_angle=None):
     """
+
     @param file: SESANS file location
     @param model: Bumps model object or model name - can be model, model_1 * model_2, and/or model_1 + model_2
@@ -55 +57 @@
             dy = err_data
             sample = Sample()
+            acceptance_angle = acceptance_angle
+            needs_all_q = acceptance_angle is not None
         data = SESANSData1D()
 

sasmodels/core.py

@@ -176 +176 @@
     return value, weight
 
-def call_kernel(kernel, pars, cutoff=0):
+def call_kernel(kernel, pars, cutoff=0, mono=False):
     """
     Call *kernel* returned from :func:`make_kernel` with parameters *pars*.
@@ -189 +189 @@
     fixed_pars = [pars.get(name, kernel.info['defaults'][name])
                   for name in kernel.fixed_pars]
-    pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
+    if mono:
+        pd_pars = [( np.array([pars[name]]), np.array([1.0]) )
+                   for name in kernel.pd_pars]
+    else:
+        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
     return kernel(fixed_pars, pd_pars, cutoff=cutoff)
 

sasmodels/direct_model.py

@@ -69 +69 @@
 
         if self.data_type == 'sesans':
+            
             q = sesans.make_q(data.sample.zacceptance, data.Rmax)
             index = slice(None, None)
@@ -77 +78 @@
                 Iq, dIq = None, None
             #self._theory = np.zeros_like(q)
-            q_vectors = [q]
+            q_vectors = [q]            
+            q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
             if not partype['orientation'] and not partype['magnetic']:
@@ -96 +98 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = res.q_calc
+            q_mono = []
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -120 +123 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = [res.q_calc]
+            q_mono = []
         else:
             raise ValueError("Unknown data type") # never gets here
@@ -125 +129 @@
         # Remember function inputs so we can delay loading the function and
         # so we can save/restore state
-        self._kernel_inputs = [v for v in q_vectors]
+        self._kernel_inputs = q_vectors
+        self._kernel_mono_inputs = q_mono
         self._kernel = None
         self.Iq, self.dIq, self.index = Iq, dIq, index
@@ -149 +154 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init
+            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type
+            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None
 
         Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
+        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None
         if self.data_type == 'sesans':
-            result = sesans.hankel(self._data.x, self._data.lam * 1e-9,
-                                   self._data.sample.thickness / 10,
-                                   self._kernel_inputs[0], Iq_calc)
+            result = sesans.transform(self._data,
+                                   self._kernel_inputs[0], Iq_calc, 
+                                   self._kernel_mono_inputs, Iq_mono)
         else:
             result = self.resolution.apply(Iq_calc)
-        return result
+        return result        
 
 

sasmodels/sesans.py

@@ -13 +13 @@
 from numpy import pi, exp
 from scipy.special import jv as besselj
-
+#import direct_model.DataMixin as model
+        
 def make_q(q_max, Rmax):
     r"""
@@ -21 +22 @@
     q_min = dq = 0.1 * 2*pi / Rmax
     return np.arange(q_min, q_max, dq)
+    
+def make_allq(data):
+    if not data.needs_all_q:
+        return []
+    elif needs_Iqxy(data):
+        # compute qx, qy
+        Qx, Qy = np.meshgrid(qx, qy)
+        return [Qx, Qy]
+    else:
+        # else only need q
+        return [q]
 
+def transform(data, q_calc, Iq_calc, qmono, Iq_mono):
+    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.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  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 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 = 
+    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"""
@@ -44 +161 @@
     """
     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