Changeset c97724e in sasmodels


Ignore:
Timestamp:
Feb 18, 2015 8:45:23 AM (10 years ago)
Author:
pkienzle
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:
e806077
Parents:
f173599
Message:

add sesans support to bumps model

Files:
2 added
1 deleted
2 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/bumps_model.py

    rf786ff3 rc97724e  
    44import sys, os 
    55import datetime 
     6 
     7from sasmodels import sesans 
    68 
    79# CRUFT python 2.6 
     
    213215    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True) 
    214216 
     217def _plot_sesans(data, theory, view): 
     218    import matplotlib.pyplot as plt 
     219    resid = (theory - data.data)/data.err_data 
     220    plt.subplot(121) 
     221    plt.errorbar(data.SElength, data.data, yerr=data.err_data) 
     222    plt.plot(data.SElength, theory, '-', hold=True) 
     223    plt.xlabel('spin echo length (A)') 
     224    plt.ylabel('polarization') 
     225    plt.subplot(122) 
     226    plt.plot(data.SElength, resid, 'x') 
     227    plt.xlabel('spin echo length (A)') 
     228    plt.ylabel('residuals') 
    215229 
    216230def _plot_result2D(data, theory, view): 
     
    230244    plt.colorbar() 
    231245 
    232 def plot_result(data, theory, view='log'): 
    233     """ 
    234     Plot the data and residuals. 
    235     """ 
    236     if hasattr(data, 'qx_data'): 
    237         _plot_result2D(data, theory, view) 
    238     else: 
    239         _plot_result1D(data, theory, view) 
    240  
    241  
    242246class BumpsModel(object): 
    243247    """ 
     
    263267        self.model = model 
    264268        self.cutoff = cutoff 
     269        if hasattr(data, 'SElength'): 
     270            self.data_type = 'sesans' 
     271        elif hasattr(data, 'qx_data'): 
     272            self.data_type = 'Iqxy' 
     273        else: 
     274            self.data_type = 'Iq' 
    265275 
    266276        partype = model.info['partype'] 
    267277 
    268278        # interpret data 
    269         if hasattr(data, 'qx_data'): 
     279        if self.data_type == 'sesans': 
     280            q = sesans.make_q(data.q_zmax, data.Rmax) 
     281            self.index = slice(None,None) 
     282            self.iq = data.data 
     283            self.diq = data.err_data 
     284            self._theory = np.zeros_like(q) 
     285            q_vectors = [q] 
     286        elif self.data_type == 'Iqxy': 
    270287            self.index = (data.mask==0) & (~np.isnan(data.data)) 
    271288            self.iq = data.data[self.index] 
     
    276293            else: 
    277294                q_vectors = [data.qx_data, data.qy_data] 
    278         else: 
     295        elif self.data_type == 'Iq': 
    279296            self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y) 
    280297            self.iq = data.y[self.index] 
     
    282299            self._theory = np.zeros_like(data.y) 
    283300            q_vectors = [data.x] 
     301        else: 
     302            raise ValueError("Unknown data type") # never gets here 
    284303 
    285304        # Remember function inputs so we can delay loading the function and 
     
    333352            self._theory[self.index] = self._fn(pars, pd_pars, self.cutoff) 
    334353            #self._theory[:] = self._fn.eval(pars, pd_pars) 
    335             self._cache['theory'] = self._theory 
     354            if self.data_type == 'sesans': 
     355                P = sesans.hankel(self.data.SElength, self.data.wavelength, 
     356                                  self.data.thickness, self._fn_inputs[0], 
     357                                  self._theory) 
     358                self._cache['theory'] = P 
     359            else: 
     360                self._cache['theory'] = self._theory 
    336361        return self._cache['theory'] 
    337362 
     
    349374 
    350375    def plot(self, view='log'): 
    351         plot_result(self.data, self.theory(), view=view) 
     376        """ 
     377        Plot the data and residuals. 
     378        """ 
     379        data, theory = self.data, self.theory() 
     380        if self.data_type == 'Iq': 
     381            _plot_result1D(data, theory, view) 
     382        elif self.data_type == 'Iqxy': 
     383            _plot_result2D(data, theory, view) 
     384        elif self.data_type == 'sesans': 
     385            _plot_sesans(data, theory, view) 
     386        else: 
     387            raise ValueError("Unknown data type") 
     388 
     389    def simulate_data(self, noise=None): 
     390        print "noise", noise 
     391        if noise is None: 
     392            noise = self.diq[self.index] 
     393        else: 
     394            noise = 0.01*noise 
     395            self.diq[self.index] = noise 
     396        y = self.theory() 
     397        y += y*np.random.randn(*y.shape)*noise 
     398        if self.data_type == 'Iq': 
     399            self.data.y[self.index] = y 
     400        elif self.data_type == 'Iqxy': 
     401            self.data.data[self.index] = y 
     402        elif self.data_type == 'sesans': 
     403            self.data.data[self.index] = y 
     404        else: 
     405            raise ValueError("Unknown model") 
    352406 
    353407    def save(self, basename): 
  • sesansdemo.py

    r10576d1 rc97724e  
    1010 
    1111# q-range parameters 
     12 
    1213q = arange(0.0003, 1.0, 0.0003);    # [nm^-1] range wide enough for  Hankel transform 
    1314dq=(q[1]-q[0])*1e9;   # [m^-1] step size in q, needed for integration 
     
    3435 
    3536clf() 
    36 subplot(1,2,1)  # plot the SANS calculation 
     37subplot(211)  # plot the SANS calculation 
    3738plot(q,I,'k') 
    3839loglog(q,I) 
     
    5556PP=exp(th*Lambda**2/4/pi/pi*(G-G[0])); 
    5657 
    57 subplot(1,2,2) 
     58subplot(212) 
    5859plot(zz,PP,'k',label="Hankel transform") # Hankel transform 1D 
    5960xlabel('spin-echo length [nm]') 
     
    6263 
    6364# Cosine transformation of 2D scattering patern 
    64 if True: 
     65if False: 
    6566    qy,qz = meshgrid(q,q) 
    6667    qr=R*sqrt(qy**2 + qz**2); # reuse variable names Hankel transform, but now 2D 
Note: See TracChangeset for help on using the changeset viewer.