Changeset 02e70ff in sasmodels for sasmodels/sesans.py

Timestamp:

Mar 17, 2016 11:03:42 AM (8 years ago)

Author:

jhbakker

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

6ddd6e0

Parents:

c970053

Message:

Beginnings of handling acceptance angle in sesans

File:

• sasmodels/sesans.py (modified) (3 diffs)

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