Changes in / [23054a1:237b7392] in sasmodels


Ignore:
Files:
1 deleted
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/sesans.py

    r94d13f1 r7fa0e08  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import j0 
    17  
    18 class SesansTransform(object): 
    19     """ 
    20     Spin-Echo SANS transform calculator.  Similar to a resolution function, 
    21     the SesansTransform object takes I(q) for the set of *q_calc* values and 
    22     produces a transformed dataset 
    23  
    24     *SElength* (A) is the set of spin-echo lengths in the measured data. 
    25  
    26     *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 
    27     echo encoding dimension (for ToF: Q of min(R) and max(lam)). 
    28  
    29     *Rmax* (A) is the maximum size sensitivity; larger radius requires more 
    30     computation time. 
    31     """ 
    32     #: SElength from the data in the original data units; not used by transform 
    33     #: but the GUI uses it, so make sure that it is present. 
    34     q = None  # type: np.ndarray 
    35  
    36     #: q values to calculate when computing transform 
    37     q_calc = None  # type: np.ndarray 
    38  
    39     # transform arrays 
    40     _H = None  # type: np.ndarray 
    41     _H0 = None # type: np.ndarray 
    42  
    43     def __init__(self, z, SElength, zaccept, Rmax): 
    44         # type: (np.ndarray, float, float) -> None 
    45         #import logging; logging.info("creating SESANS transform") 
    46         self.q = z 
    47         self._set_hankel(SElength, zaccept, Rmax) 
    48  
    49     def apply(self, Iq): 
    50         # tye: (np.ndarray) -> np.ndarray 
    51         G0 = np.dot(self._H0, Iq) 
    52         G = np.dot(self._H.T, Iq) 
    53         P = G - G0 
    54         return P 
    55  
    56     def _set_hankel(self, SElength, zaccept, Rmax): 
    57         # type: (np.ndarray, float, float) -> None 
    58         # Force float32 arrays, otherwise run into memory problems on some machines 
    59         SElength = np.asarray(SElength, dtype='float32') 
    60  
    61         #Rmax = #value in text box somewhere in FitPage? 
    62         q_max = 2*pi / (SElength[1] - SElength[0]) 
    63         q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 
    64         q = np.arange(q_min, q_max, q_min, dtype='float32') 
    65         dq = q_min 
    66  
    67         H0 = np.float32(dq/(2*pi)) * q 
    68  
    69         repq = np.tile(q, (SElength.size, 1)).T 
    70         repSE = np.tile(SElength, (q.size, 1)) 
    71         H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 
    72  
    73         self.q_calc = q 
    74         self._H, self._H0 = H, H0 
     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 TracChangeset for help on using the changeset viewer.