Changeset 02e70ff in sasmodels


Ignore:
Timestamp:
Mar 17, 2016 1:03:42 PM (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

Files:
4 edited

Legend:

Unmodified
Added
Removed

  • example/sesansfit.py

    r84db7a5 r02e70ff  
    88    HAS_CONVERTER = False 
    99 
     10 
    1011def get_bumps_model(model_name): 
    1112    kernel = core.load_model(model_name) 
    1213    model = bumps_model.Model(kernel) 
    1314    return model 
     15 
     16def sesans_fit(file, model_name, initial_vals={}, custom_params={}, param_range=[], acceptance_angle=None): 
     17    """ 
    1418 
    1519def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[]): 
    … …  
    5559            dy = err_data 
    5660            sample = Sample() 
     61            acceptance_angle = acceptance_angle 
     62            needs_all_q = acceptance_angle is not None 
    5763        data = SESANSData1D() 
    5864 

  • sasmodels/core.py

    r7b3e62c r02e70ff  
    176176    return value, weight 
    177177 
    178 def call_kernel(kernel, pars, cutoff=0): 
     178def call_kernel(kernel, pars, cutoff=0, mono=False): 
    179179    """ 
    180180    Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
    … …  
    189189    fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    190190                  for name in kernel.fixed_pars] 
    191     pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
     191    if mono: 
     192        pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
     193                   for name in kernel.pd_pars] 
     194    else: 
     195        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
    192196    return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
    193197 

  • sasmodels/direct_model.py

    r17bbadd r02e70ff  
    6969 
    7070        if self.data_type == 'sesans': 
     71             
    7172            q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
    7273            index = slice(None, None) 
    … …  
    7778                Iq, dIq = None, None 
    7879            #self._theory = np.zeros_like(q) 
    79             q_vectors = [q] 
     80            q_vectors = [q]             
     81            q_mono = sesans.make_all_q(data) 
    8082        elif self.data_type == 'Iqxy': 
    8183            if not partype['orientation'] and not partype['magnetic']: 
    … …  
    9698            #self._theory = np.zeros_like(self.Iq) 
    9799            q_vectors = res.q_calc 
     100            q_mono = [] 
    98101        elif self.data_type == 'Iq': 
    99102            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    … …  
    120123            #self._theory = np.zeros_like(self.Iq) 
    121124            q_vectors = [res.q_calc] 
     125            q_mono = [] 
    122126        else: 
    123127            raise ValueError("Unknown data type") # never gets here 
    … …  
    125129        # Remember function inputs so we can delay loading the function and 
    126130        # so we can save/restore state 
    127         self._kernel_inputs = [v for v in q_vectors] 
     131        self._kernel_inputs = q_vectors 
     132        self._kernel_mono_inputs = q_mono 
    128133        self._kernel = None 
    129134        self.Iq, self.dIq, self.index = Iq, dIq, index 
    … …  
    149154    def _calc_theory(self, pars, cutoff=0.0): 
    150155        if self._kernel is None: 
    151             self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init 
     156            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
     157            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None 
    152158 
    153159        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
     160        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None 
    154161        if self.data_type == 'sesans': 
    155             result = sesans.hankel(self._data.x, self._data.lam * 1e-9, 
    156                                    self._data.sample.thickness / 10, 
    157                                    self._kernel_inputs[0], Iq_calc) 
     162            result = sesans.transform(self._data, 
     163                                   self._kernel_inputs[0], Iq_calc,  
     164                                   self._kernel_mono_inputs, Iq_mono) 
    158165        else: 
    159166            result = self.resolution.apply(Iq_calc) 
    160         return result 
     167        return result         
    161168 
    162169 

  • sasmodels/sesans.py

    r190fc2b r02e70ff  
    1313from numpy import pi, exp 
    1414from scipy.special import jv as besselj 
    15  
     15#import direct_model.DataMixin as model 
     16         
    1617def make_q(q_max, Rmax): 
    1718    r""" 
    … …  
    2122    q_min = dq = 0.1 * 2*pi / Rmax 
    2223    return np.arange(q_min, q_max, dq) 
     24     
     25def make_allq(data): 
     26    if not data.needs_all_q: 
     27        return [] 
     28    elif needs_Iqxy(data): 
     29        # compute qx, qy 
     30        Qx, Qy = np.meshgrid(qx, qy) 
     31        return [Qx, Qy] 
     32    else: 
     33        # else only need q 
     34        return [q] 
    2335 
     36def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
     37    nqmono = len(qmono) 
     38    if nqmono == 0: 
     39        result = call_hankel(data, q_calc, Iq_calc) 
     40    elif nqmono == 1: 
     41        q = qmono[0] 
     42        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
     43    else: 
     44        Qx, Qy = [qmono[0], qmono[1]] 
     45        Qx = np.reshape(Qx, nqx, nqy) 
     46        Qy = np.reshape(Qy, nqx, nqy) 
     47        Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
     48        qx = Qx[0, :] 
     49        qy = Qy[:, 0] 
     50        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     51 
     52    return result 
     53 
     54def call_hankel(data, q_calc, Iq_calc): 
     55    return hankel(data.x, data.lam * 1e-9, 
     56                  data.sample.thickness / 10, 
     57                  q_calc, Iq_calc) 
     58   
     59def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
     60    return hankel(data.x, data.lam * 1e-9, 
     61                  data.sample.thickness / 10, 
     62                  q_calc, Iq_calc) 
     63                   
     64def Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
     65    return hankel(data.x, data.y data.lam * 1e-9, 
     66                  data.sample.thickness / 10, 
     67                  q_calc, Iq_calc) 
     68                         
     69def TotalScatter(model, parameters):  #Work in progress!! 
     70#    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) 
     71    allq = np.linspace(0,4*pi/wavelength,1000) 
     72    allIq =  
     73    integral = allq*allIq 
     74     
     75 
     76 
     77def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
     78#============================================================================== 
     79#     2D Cosine Transform if "wavelength" is a vector 
     80#============================================================================== 
     81#allq is the q-space needed to create the total scattering cross-section 
     82 
     83    Gprime = np.zeros_like(wavelength, 'd') 
     84    s = np.zeros_like(wavelength, 'd') 
     85    sd = np.zeros_like(wavelength, 'd') 
     86    Gprime = np.zeros_like(wavelength, 'd') 
     87    f = np.zeros_like(wavelength, 'd') 
     88       for i, wavelength_i in enumerate(wavelength): 
     89            z = magfield*wavelength_i 
     90            allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
     91            allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
     92            alldq = (allq[1]-allq[0])*1e10 
     93            sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
     94            s[i]=1-exp(-sigma) 
     95            for j, Iqy_j, qy_j in enumerate(qy): 
     96                for k, Iqz_k, qz_k in enumerate(qz): 
     97                    Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
     98                    q = np.sqrt(qy_j^2 + qz_k^2) 
     99                    Gintegral = Iq*cos(z*Qz_k) 
     100                    Gprime[i] += Gintegral 
     101    #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
     102    #                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
     103    #                For now, work with standard 2-phase scatter 
     104                    
     105                     
     106                    sd[i] += Iq 
     107            f[i] = 1-s[i]+sd[i] 
     108            P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]         
     109 
     110 
     111 
     112 
     113def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
     114#============================================================================== 
     115#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
     116#============================================================================== 
     117#acceptq is the q-space needed to create limited acceptance effect 
     118    SElength= wavelength*magfield 
     119    G = np.zeros_like(SElength, 'd') 
     120    threshold=2*pi*theta/wavelength 
     121        for i, SElength_i in enumerate(SElength): 
     122            allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
     123            allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
     124            alldq = (allq[1]-allq[0])*1e10 
     125            sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
     126            s[i]=1-exp(-sigma) 
     127             
     128            dq = (q[1]-q[0])*1e10 
     129            a = (x<threshold) 
     130            acceptq = a*q 
     131            acceptIq = a*Iq 
     132        
     133            G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
     134                 
     135    #        G[i]=np.sum(integral) 
     136         
     137        G *= dq*1e10*2*pi 
     138     
     139        P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
     140     
    24141def hankel(SElength, wavelength, thickness, q, Iq): 
    25142    r""" 
    … …  
    44161    """ 
    45162    G = np.zeros_like(SElength, 'd') 
     163#============================================================================== 
     164#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
     165#============================================================================== 
    46166    for i, SElength_i in enumerate(SElength): 
    47167        integral = besselj(0, q*SElength_i)*Iq*q 
Note: See TracChangeset for help on using the changeset viewer.