source: sasmodels/sasmodels/sesans.py @ b397165

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b397165 was b397165, checked in by GitHub <noreply@…>, 7 years ago

Revert "Jurtest2"

  • Property mode set to 100644
File size: 8.4 KB
Line 
1"""
2Conversion of scattering cross section from SANS (I(q), or rather, ds/dO) in absolute
3units (cm-1)into SESANS correlation function G using a Hankel transformation, then converting
4the SESANS correlation function into polarisation from the SESANS experiment
5
6Everything is in units of metres except specified otherwise (NOT TRUE!!!)
7Everything is in conventional units (nm for spin echo length)
8
9Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013
10"""
11
12from __future__ import division
13
14import numpy as np  # type: ignore
15from numpy import pi, exp  # type: ignore
16from scipy.special import jv as besselj
17#import direct_model.DataMixin as model
18       
19def make_q(q_max, Rmax):
20    r"""
21    Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$
22    to $q_max$. This is the integration range of the Hankel transform; bigger range and
23    more points makes a better numerical integration.
24    Smaller q_min will increase reliable spin echo length range.
25    Rmax is the "radius" of the largest expected object and can be set elsewhere.
26    q_max is determined by the acceptance angle of the SESANS instrument.
27    """
28    from sas.sascalc.data_util.nxsunit import Converter
29
30    q_min = dq = 0.1 * 2*pi / Rmax
31    return np.arange(q_min,
32                     Converter(q_max[1])(q_max[0],
33                                         units="1/A"),
34                     dq)
35   
36def make_all_q(data):
37    """
38    Return a $q$ vector suitable for calculating the total scattering cross section for
39    calculating the effect of finite acceptance angles on Time of Flight SESANS instruments.
40    If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed.
41    If the instrument has a rectangular acceptance, 2 all_q vectors are needed.
42    If the instrument has a circular acceptance, 1 all_q vector is needed
43   
44    """
45    if not data.has_no_finite_acceptance:
46        return []
47    elif data.has_yz_acceptance(data):
48        # compute qx, qy
49        Qx, Qy = np.meshgrid(qx, qy)
50        return [Qx, Qy]
51    else:
52        # else only need q
53        # data.has_z_acceptance
54        return [q]
55
56def transform(data, q_calc, Iq_calc, qmono, Iq_mono):
57    """
58    Decides which transform type is to be used, based on the experiment data file contents (header)
59    (2016-03-19: currently controlled from parameters script)
60    nqmono is the number of q vectors to be used for the detector integration
61    """
62    nqmono = len(qmono)
63    if nqmono == 0:
64        result = call_hankel(data, q_calc, Iq_calc)
65    elif nqmono == 1:
66        q = qmono[0]
67        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono)
68    else:
69        Qx, Qy = [qmono[0], qmono[1]]
70        Qx = np.reshape(Qx, nqx, nqy)
71        Qy = np.reshape(Qy, nqx, nqy)
72        Iq_mono = np.reshape(Iq_mono, nqx, nqy)
73        qx = Qx[0, :]
74        qy = Qy[:, 0]
75        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono)
76
77    return result
78
79def call_hankel(data, q_calc, Iq_calc):
80    return hankel((data.x, data.x_unit),
81                  (data.lam, data.lam_unit),
82                  (data.sample.thickness,
83                   data.sample.thickness_unit),
84                  q_calc, Iq_calc)
85 
86def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono):
87    return hankel(data.x, data.lam * 1e-9,
88                  data.sample.thickness / 10,
89                  q_calc, Iq_calc)
90                 
91def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono):
92    return hankel(data.x, data.y, data.lam * 1e-9,
93                  data.sample.thickness / 10,
94                  q_calc, Iq_calc)
95                       
96def TotalScatter(model, parameters):  #Work in progress!!
97#    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)
98    allq = np.linspace(0,4*pi/wavelength,1000)
99    allIq = 1
100    integral = allq*allIq
101   
102
103
104def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still
105#==============================================================================
106#     2D Cosine Transform if "wavelength" is a vector
107#==============================================================================
108#allq is the q-space needed to create the total scattering cross-section
109
110    Gprime = np.zeros_like(wavelength, 'd')
111    s = np.zeros_like(wavelength, 'd')
112    sd = np.zeros_like(wavelength, 'd')
113    Gprime = np.zeros_like(wavelength, 'd')
114    f = np.zeros_like(wavelength, 'd')
115    for i, wavelength_i in enumerate(wavelength):
116        z = magfield*wavelength_i
117        allq=np.linspace() #for calculating the Q-range of the  scattering power integral
118        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
119        alldq = (allq[1]-allq[0])*1e10
120        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
121        s[i]=1-exp(-sigma)
122        for j, Iqy_j, qy_j in enumerate(qy):
123            for k, Iqz_k, qz_k in enumerate(qz):
124                Iq = np.sqrt(Iqy_j^2+Iqz_k^2)
125                q = np.sqrt(qy_j^2 + qz_k^2)
126                Gintegral = Iq*cos(z*Qz_k)
127                Gprime[i] += Gintegral
128#                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i]
129#                s[i] += 1-exp(Totalscatter(modelname)*thickness)
130#                For now, work with standard 2-phase scatter
131
132
133                sd[i] += Iq
134        f[i] = 1-s[i]+sd[i]
135        P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]
136
137
138
139
140def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname):
141#==============================================================================
142#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS
143#==============================================================================
144#acceptq is the q-space needed to create limited acceptance effect
145    SElength= wavelength*magfield
146    G = np.zeros_like(SElength, 'd')
147    threshold=2*pi*theta/wavelength
148    for i, SElength_i in enumerate(SElength):
149        allq=np.linspace() #for calculating the Q-range of the  scattering power integral
150        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
151        alldq = (allq[1]-allq[0])*1e10
152        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
153        s[i]=1-exp(-sigma)
154
155        dq = (q[1]-q[0])*1e10
156        a = (x<threshold)
157        acceptq = a*q
158        acceptIq = a*Iq
159
160        G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq)
161
162#        G[i]=np.sum(integral)
163
164    G *= dq*1e10*2*pi
165
166    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
167   
168def hankel(SElength, wavelength, thickness, q, Iq):
169    r"""
170    Compute the expected SESANS polarization for a given SANS pattern.
171
172    Uses the hankel transform followed by the exponential.  The values for *zz*
173    (or spin echo length, or delta), wavelength and sample thickness should
174    come from the dataset.  $q$ should be chosen such that the oscillations
175    in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$).
176
177    *SElength* [A] is the set of $z$ points at which to compute the
178    Hankel transform
179
180    *wavelength* [m]  is the wavelength of each individual point *zz*
181
182    *thickness* [cm] is the sample thickness.
183
184    *q* [A$^{-1}$] is the set of $q$ points at which the model has been
185    computed. These should be equally spaced.
186
187    *I* [cm$^{-1}$] is the value of the SANS model at *q*
188    """
189
190    from sas.sascalc.data_util.nxsunit import Converter
191    wavelength = Converter(wavelength[1])(wavelength[0],"A")
192    thickness = Converter(thickness[1])(thickness[0],"A")
193    Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters
194    SElength = Converter(SElength[1])(SElength[0],"A")
195
196    G = np.zeros_like(SElength, 'd')
197#==============================================================================
198#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS
199#==============================================================================
200    for i, SElength_i in enumerate(SElength):
201        integral = besselj(0, q*SElength_i)*Iq*q
202        G[i] = np.sum(integral)
203    G0 = np.sum(Iq*q)
204
205    # [m^-1] step size in q, needed for integration
206    dq = (q[1]-q[0])
207
208    # integration step, convert q into [m**-1] and 2 pi circle integration
209    G *= dq*2*pi
210    G0 = np.sum(Iq*q)*dq*2*np.pi
211
212    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0))
213
214    return P
Note: See TracBrowser for help on using the repository browser.