Changeset 54f1d96 in sasmodels


Ignore:
Timestamp:
Jan 17, 2017 2:50:42 PM (6 years ago)
Author:
jhbakker
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
2f8ee5f
Parents:
daeef4c
Message:

Manually added in all SESANS modifications from Jurtest branch because of
corrupted git history…

Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/sesans.py

    rb397165 r54f1d96  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import jv as besselj 
    17 #import direct_model.DataMixin as model 
     16from scipy.special import j0 
     17#from mpmath import j0 as j0 
    1818         
    19 def 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 
     19class SesansTransform(object): 
     20    #: Set of spin-echo lengths in the measured data 
     21    SE = None  # type: np.ndarray 
     22    #: Maximum acceptance of scattering vector in the spin echo encoding dimension (for ToF: Q of min(R) and max(lam)) 
     23    zaccept = None # type: float 
     24    #: Maximum size sensitivity; larger radius requires more computation 
     25    Rmax = None  # type: float 
     26    #: q values to calculate when computing transform 
     27    q = None  # type: np.ndarray 
    2928 
    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      
    36 def 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] 
     29    # transform arrays 
     30    _H = None  # type: np.ndarray 
     31    _H0 = None # type: np.ndarray 
    5532 
    56 def 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) 
     33    def set_transform(self, SE, zaccept, Rmax): 
     34        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: 
     35            self.SE, self.zaccept, self.Rmax = SE, zaccept, Rmax 
     36            self._set_q() 
     37            self._set_hankel() 
    7638 
    77     return result 
     39    def apply(self, Iq): 
     40        G0 = np.dot(self._H0, Iq) 
     41        G = np.dot(self._H.T, Iq) 
     42        P = G - G0 
     43        return P 
    7844 
    79 def 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    
    86 def 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                    
    91 def 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                          
    96 def 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      
     45    def _set_q(self): 
     46        #q_min = dq = 0.1 * 2*pi / self.Rmax 
    10247 
     48        q_max = 2*pi / (self.SE[1]-self.SE[0]) 
     49        q_min = dq = 0.1 *2*pi / (np.size(self.SE) * self.SE[-1]) 
    10350 
    104 def 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 
     51        #q_min = dq = q_max / 100000 
     52        q=np.arange(q_min, q_max, q_min) 
     53        self.q = q 
     54        self.dq = dq 
    10955 
    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 
     56    def _set_hankel(self): 
     57        #Rmax = #value in text box somewhere in FitPage? 
     58        q = self.q 
     59        dq = self.dq 
     60        SElength = self.SE 
    13161 
     62        H0 = dq / (2 * pi) * q 
     63        q=np.array(q,dtype='float32') 
     64        SElength=np.array(SElength,dtype='float32') 
    13265 
    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] 
     66        # Using numpy tile, dtype is conserved 
     67        repq=np.tile(q,(SElength.size,1)) 
     68        repSE=np.tile(SElength,(q.size,1)) 
     69        H = dq / (2 * pi) * j0(repSE*repq.T)*repq.T 
    13670 
     71        # Using numpy meshgrid - meshgrid produces float64 from float32 inputs! Problem for 32-bit OS: Memerrors! 
     72        #H0 = dq / (2 * pi) * q 
     73        #repSE, repq = np.meshgrid(SElength, q) 
     74        #repq=np.array(repq,dtype='float32') 
     75        #repSE=np.array(repSE,dtype='float32') 
     76        #H = dq / (2 * pi) * j0(repSE*repq)*repq 
    13777 
     78        self._H, self._H0 = H, H0 
    13879 
    139  
    140 def 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      
    168 def 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 
     80class SESANS1D(SesansTransform): 
     81    def __init__(self, data, _H0, _H, q_calc): 
     82        # x values of the data (Sasview requires these to be named "q") 
     83        self.q = data.x 
     84        self._H0 = _H0 
     85        self._H = _H 
     86        # Pysmear does some checks on the smearer object, these checks require the "data" object... 
     87        self.data=data 
     88        # q values of the SAS model 
     89        self.q_calc = q_calc # These are the MODEL's q values used by the smearer (in this case: the Hankel transform) 
     90    def apply(self, theory): 
     91        return SesansTransform.apply(self,theory) 
Note: See TracChangeset for help on using the changeset viewer.