Changeset eec5fd6 in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 4:58:02 AM (9 years ago)
Author:
richardh
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:
ead4bd5
Parents:
4f4e0d5 (diff), 4554131 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed

  • doc/genmodel.py

    raa2edb2 rc094758  
    1010# Convert ../sasmodels/models/name.py to name 
    1111model_name = os.path.basename(sys.argv[1])[:-3] 
     12model_info = core.load_model_info(model_name) 
     13model = core.build_model(model_info) 
    1214 
    1315# Load the doc string from the module definition file and store it in rst 
    14 docstr = generate.make_doc(core.load_model_info(model_name)) 
    15      
    16 # Generate automatically plot of the model and add it to rst documentation 
     16docstr = generate.make_doc(model_info) 
    1717 
    18 info = core.load_model_info(model_name) 
    1918 
    2019# Calculate 1D curve for default parameters 
    21 pars = dict((p[0], p[2]) for p in info['parameters']) 
     20pars = dict((p.name, p.default) for p in model_info['parameters']) 
    2221 
    2322# Plotting ranges and options 
    2423opts = { 
    25         'xscale'    : 'log', 
    26         'yscale'    : 'log' if not info['structure_factor'] else 'linear', 
    27         'qmin'      : 0.001, 
    28         'qmax'      : 1.0, 
    29         'nq'        : 1000, 
    30         'nq2d'      : 100, 
     24    'xscale'    : 'log', 
     25    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     26    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     27    'q_min'     : 0.001, 
     28    'q_max'     : 1.0, 
     29    'nq'        : 1000, 
     30    'nq2d'      : 400, 
     31    'vmin'      : 1e-3,  # floor for the 2D data results 
     32    'qx_max'    : 0.5, 
    3133} 
    3234 
    33 qmin, qmax, nq = opts['qmin'], opts['qmax'], opts['nq'] 
    34 qmin = math.log10(qmin) 
    35 qmax = math.log10(qmax) 
    36 q = np.logspace(qmin, qmax, nq) 
    37 data = empty_data1D(q) 
    38 model = core.load_model(model_name) 
    39 calculator = DirectModel(data, model) 
    40 Iq1D = calculator() 
    4135 
    42 # TO DO: Generation of 2D plots  
    43 # Problem in sasmodels.direct_model._calc_theory 
    44 # There self._kernel.q_input.nq  gets a value of 0 in the 2D case 
    45 # and returns a 0 numpy array (it does not call the C code) 
     36def plot_1d(model, opts, ax): 
     37    q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq'] 
     38    q_min = math.log10(q_min) 
     39    q_max = math.log10(q_max) 
     40    q = np.logspace(q_min, q_max, nq) 
     41    data = empty_data1D(q) 
     42    calculator = DirectModel(data, model) 
     43    Iq1D = calculator() 
    4644 
    47 # If 2D model, compute 2D image 
    48 #if info['has_2d'] != []: 
    49 #    qmax, nq2d = opts['qmax'], opts['nq2d'] 
    50 #    data2d = empty_data2D(np.linspace(-qmax, qmax, nq2d), resolution=0.0)  
    51 #    #model = core.load_model(model_name) 
    52 #    calculator = DirectModel(data2d, model) 
    53 #    Iq2D = calculator() 
     45    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name']) 
     46    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
     47    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
     48    ax.set_xscale(opts['xscale']) 
     49    ax.set_yscale(opts['yscale']) 
     50    #ax.legend(loc='best') 
     51 
     52def plot_2d(model, opts, ax): 
     53    qx_max, nq2d = opts['qx_max'], opts['nq2d'] 
     54    q = np.linspace(-qx_max, qx_max, nq2d) 
     55    data2d = empty_data2D(q, resolution=0.0) 
     56    calculator = DirectModel(data2d, model) 
     57    Iq2D = calculator() #background=0) 
     58    Iq2D = Iq2D.reshape(nq2d, nq2d) 
     59    if opts['zscale'] == 'log': 
     60        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 
     61    h = ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='upper', 
     62           extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=ice_cm()) 
     63    # , vmin=vmin, vmax=vmax) 
     64    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 
     65    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 
     66 
     67def ice_cm(): 
     68    from matplotlib._cm import _Blues_data 
     69    from matplotlib import colors 
     70    from matplotlib import rcParams 
     71    def from_white(segments): 
     72        scale = 1.0/segments[0][1] 
     73        return [(k, v*scale, w*scale) for k, v, w in segments] 
     74    ice_data = dict((k,from_white(v)) for k,v in _Blues_data.items()) 
     75    ice = colors.LinearSegmentedColormap("ice", ice_data, rcParams['image.lut']) 
     76    return ice 
     77 
    5478 
    5579# Generate image (comment IF for 1D/2D for the moment) and generate only 1D 
    56 #if info['has_2d'] == []: 
    57 #    fig = plt.figure() 
    58 #    ax = fig.add_subplot(1,1,1) 
    59 #    ax.plot(q, Iq1D, color='blue', lw=2, label=model_name) 
    60 #    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
    61 #    ax.set_xscale(opts['xscale'])    
    62 #    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
    63 #    ax.set_yscale(opts['yscale'])   
    64 #    ax.legend() 
    65 #else: 
    66 #    # need figure with 1D + 2D 
    67 #    pass 
    68 fig = plt.figure() 
    69 ax = fig.add_subplot(1,1,1) 
    70 ax.plot(q, Iq1D, color='blue', lw=2, label=info['name']) 
    71 ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
    72 ax.set_xscale(opts['xscale'])    
    73 ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
    74 ax.set_yscale(opts['yscale'])   
    75 ax.legend() 
    76   
     80fig_height = 3.0 # in 
     81fig_left = 0.6 # in 
     82fig_right = 0.5 # in 
     83fig_top = 0.6*0.25 # in 
     84fig_bottom = 0.6*0.75 
     85if model_info['has_2d']: 
     86    plot_height = fig_height - (fig_top+fig_bottom) 
     87    plot_width = plot_height 
     88    fig_width = 2*(plot_width + fig_left + fig_right) 
     89    aspect = (fig_width, fig_height) 
     90    ratio = aspect[0]/aspect[1] 
     91    ax_left = fig_left/fig_width 
     92    ax_bottom = fig_bottom/fig_height 
     93    ax_height = plot_height/fig_height 
     94    ax_width = ax_height/ratio # square axes 
     95    fig = plt.figure(figsize=aspect) 
     96    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 
     97    plot_2d(model, opts, ax2d) 
     98    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     99    #ax.set_aspect('square') 
     100else: 
     101    plot_height = fig_height - (fig_top+fig_bottom) 
     102    plot_width = (1+np.sqrt(5))/2*fig_height 
     103    fig_width = plot_width + fig_left + fig_right 
     104    ax_left = fig_left/fig_width 
     105    ax_bottom = fig_bottom/fig_height 
     106    ax_width = plot_width/fig_width 
     107    ax_height = plot_height/fig_height 
     108    aspect = (fig_width, fig_height) 
     109    fig = plt.figure(figsize=aspect) 
     110    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     111plot_1d(model, opts, ax1d) 
    77112 
    78113# Save image in model/img 
    79114figname = model_name + '_autogenfig.png' 
    80115filename = os.path.join('model', 'img', figname) 
    81 plt.savefig(filename) 
     116plt.savefig(filename, bbox_inches='tight') 
     117#print "figure saved in",filename 
    82118 
    83119# Auto caption for figure 
    84120captionstr = '\n' 
    85 captionstr += '.. figure:: img/' + model_name + '_autogenfig.png\n' 
     121captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n' 
    86122captionstr += '\n' 
    87 #if info['has_2d'] == []: 
    88 #    captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    89 #else: 
    90 #    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n' 
    91123captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    92124captionstr += '\n' 
    … …  
    102134else: 
    103135    print '------------------------------------------------------------------' 
    104     print 'References NOT FOUND for model: ', model_name 
     136    print 'References NOT FOUND for model: ', model_info['id'] 
    105137    print '------------------------------------------------------------------' 
    106138    docstr = docstr + captionstr 

  • example/sesansfit.py

    r84db7a5 r4554131  
    11from bumps.names import * 
    2 from sasmodels import core, bumps_model, sesans 
     2from sas import core, bumps_model, sesans 
    33 
    44HAS_CONVERTER = True 
    … …  
    88    HAS_CONVERTER = False 
    99 
     10 
    1011def get_bumps_model(model_name): 
    1112    kernel = core.load_model(model_name) 
    … …  
    1314    return model 
    1415 
    15 def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[]): 
     16def sesans_fit(file, model_name, initial_vals={}, custom_params={}, param_range=[], acceptance_angle=None): 
    1617    """ 
     18 
    1719    @param file: SESANS file location 
    1820    @param model: Bumps model object or model name - can be model, model_1 * model_2, and/or model_1 + model_2 
    … …  
    5557            dy = err_data 
    5658            sample = Sample() 
     59            acceptance_angle = acceptance_angle 
     60            needs_all_q = acceptance_angle is not None 
    5761        data = SESANSData1D() 
    5862 

  • explore/J1c.py

    r0a6da3c rcbd37a7  
    99 
    1010 
    11 SHOW_DIFF = True  # True if show diff rather than function value 
     11SHOW_DIFF = True # True if show diff rather than function value 
     12#SHOW_DIFF = False # True if show diff rather than function value 
    1213LINEAR_X = False  # True if q is linearly spaced instead of log spaced 
     14#LINEAR_X = True # True if q is linearly spaced instead of log spaced 
     15FUNCTION = "2*J1(x)/x" 
    1316 
    14 def mp_J1c(vec, bits=500): 
     17def mp_fn(vec, bits=500): 
    1518    """ 
    1619    Direct calculation using sympy multiprecision library. 
    1720    """ 
    1821    with mp.workprec(bits): 
    19         return [_mp_J1c(mp.mpf(x)) for x in vec] 
     22        return [_mp_fn(mp.mpf(x)) for x in vec] 
    2023 
    21 def _mp_J1c(x): 
     24def _mp_fn(x): 
    2225    """ 
    23     Helper funciton for mp_j1c 
     26    Actual function that gets evaluated.  The caller just vectorizes. 
    2427    """ 
    2528    return mp.mpf(2)*mp.j1(x)/x 
    2629 
    27 def np_J1c(x, dtype): 
     30def np_fn(x, dtype): 
    2831    """ 
    2932    Direct calculation using scipy. 
    … …  
    3336    return np.asarray(2, dtype)*J1(x)/x 
    3437 
    35 def cephes_J1c(x, dtype, n): 
     38def sasmodels_fn(x, dtype, platform='ocl'): 
    3639    """ 
    3740    Calculation using pade approximant. 
    3841    """ 
    39     f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    40     x = np.asarray(x, dtype) 
    41     ans = np.empty_like(x) 
    42     ax = abs(x) 
    43  
    44     # Branch a 
    45     a_idx = ax < f(8.0) 
    46     a_xsq = x[a_idx]**2 
    47     a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606))) 
    48     a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0))) 
    49     a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq) 
    50     a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq) 
    51     ans[a_idx] = f(2.0)*a_ans1/a_ans2 
    52  
    53     # Branch b 
    54     b_idx = ~a_idx 
    55     b_ax = ax[b_idx] 
    56     b_x = x[b_idx] 
    57  
    58     b_y = f(64.0)/(b_ax**2) 
    59     b_xx = b_ax - f(2.356194491) 
    60     b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6))) 
    61     b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6))) 
    62     b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y) 
    63     b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y) 
    64     b_sn, b_cn = np.sin(b_xx), np.cos(b_xx) 
    65     ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1 - (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x 
    66  
    67     return ans 
    68  
    69 def div_J1c(x, dtype): 
    70     f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    71     x = np.asarray(x, dtype) 
    72     return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype) 
    73  
    74 def _J1(x, f): 
    75     ax = abs(x) 
    76     if ax < f(8.0): 
    77         y = x*x 
    78         ans1 = x*(f(72362614232.0) 
    79                   + y*(f(-7895059235.0) 
    80                   + y*(f(242396853.1) 
    81                   + y*(f(-2972611.439) 
    82                   + y*(f(15704.48260) 
    83                   + y*(f(-30.16036606))))))) 
    84         ans2 = (f(144725228442.0) 
    85                   + y*(f(2300535178.0) 
    86                   + y*(f(18583304.74) 
    87                   + y*(f(99447.43394) 
    88                   + y*(f(376.9991397) 
    89                   + y))))) 
    90         return ans1/ans2 
    91     else: 
    92         y = f(64.0)/(ax*ax) 
    93         xx = ax - f(2.356194491) 
    94         ans1 = (f(1.0) 
    95                   + y*(f(0.183105e-2) 
    96                   + y*(f(-0.3516396496e-4) 
    97                   + y*(f(0.2457520174e-5) 
    98                   + y*f(-0.240337019e-6))))) 
    99         ans2 = (f(0.04687499995) 
    100                   + y*(f(-0.2002690873e-3) 
    101                   + y*(f(0.8449199096e-5) 
    102                   + y*(f(-0.88228987e-6) 
    103                   + y*f(0.105787412e-6))))) 
    104         sn, cn = np.sin(xx), np.cos(xx) 
    105         ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) 
    106         return -ans if (x < f(0.0)) else ans 
     42    from sasmodels import core, data, direct_model 
     43    model = core.load_model('bessel', dtype=dtype) 
     44    calculator = direct_model.DirectModel(data.empty_data1D(x), model) 
     45    return calculator(background=0) 
    10746 
    10847def plotdiff(x, target, actual, label): 
    … …  
    11352    """ 
    11453    if SHOW_DIFF: 
    115         err = np.clip(abs((target-actual)/target), 0, 1) 
     54        err = abs((target-actual)/target) 
     55        #err = np.clip(err, 0, 1) 
    11656        pylab.loglog(x, err, '-', label=label) 
    11757    else: 
    … …  
    11959        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label) 
    12060 
    121 def compare(x, precision): 
     61def compare(x, precision, target): 
    12262    r""" 
    12363    Compare the different computation methods using the given precision. 
    12464    """ 
    125     target = np.asarray(mp_J1c(x), 'double') 
    126     #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits') 
    127     plotdiff(x, target, np_J1c(x, precision), 'direct '+precision) 
    128     plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision) 
    129     #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision) 
    130     #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision) 
     65    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits') 
     66    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision) 
     67    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision) 
    13168    pylab.xlabel("qr (1/Ang)") 
    13269    if SHOW_DIFF: 
    13370        pylab.ylabel("relative error") 
    13471    else: 
    135         pylab.ylabel("2 J1(x)/x") 
     72        pylab.ylabel(FUNCTION) 
    13673        pylab.semilogx(x, target,  '-', label="true value") 
    13774    if LINEAR_X: 
    … …  
    14784    else: 
    14885        qr = np.logspace(-3,5,400) 
     86    target = np.asarray(mp_fn(qr), 'double') 
    14987    pylab.subplot(121) 
    150     compare(qr, 'single') 
     88    compare(qr, 'single', target) 
    15189    pylab.legend(loc='best') 
    15290    pylab.subplot(122) 
    153     compare(qr, 'double') 
     91    compare(qr, 'double', target) 
    15492    pylab.legend(loc='best') 
    155     pylab.suptitle('2 J1(x)/x') 
     93    pylab.suptitle(FUNCTION) 
    15694 
    15795if __name__ == "__main__": 

  • 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/data.py

    r84db7a5 rc094758  
    178178        self.qy_data, self.dqy_data = y, dy 
    179179        self.data, self.err_data = z, dz 
    180         self.mask = (~np.isnan(z) if z is not None 
    181                      else np.ones_like(x) if x is not None 
     180        self.mask = (np.isnan(z) if z is not None 
     181                     else np.zeros_like(x, dtype='bool') if x is not None 
    182182                     else None) 
    183183        self.q_data = np.sqrt(x**2 + y**2) 

  • 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/kernelcl.py

    r17bbadd rc094758  
    367367        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
    368368        context = env.get_context(self.dtype) 
     369        self.global_size = [self.q_vectors[0].size] 
     370        #print("creating inputs of size", self.global_size) 
    369371        self.q_buffers = [ 
    370372            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
    371373            for q in self.q_vectors 
    372374        ] 
    373         self.global_size = [self.q_vectors[0].size] 
    374375 
    375376    def release(self): 

  • sasmodels/models/bessel.py

    r07142f3 rcbd37a7  
    6767#Bessel 
    6868parameters = [ 
     69    ["ignored", "", 0.0, [-inf, inf], "", "no parameterless functions"], 
    6970             ] 
    7071 
    71 source = ["lib/polevl.c", "lib/j1d.c"] 
     72source = ["lib/polevl.c", "lib/j1_cephes.c"] 
    7273 
    7374# No volume normalization despite having a volume parameter 
    … …  
    7778 
    7879Iq = """ 
    79     return j1(q); 
     80    return 2.0*j1(q)/q; 
    8081    """ 
    8182 

  • sasmodels/models/lib/j0_cephes.c

    rbfef528 r094e320  
    4444 */ 
    4545 
    46 /*                                                      y0.c 
    47  * 
    48  *      Bessel function of the second kind, order zero 
    49  * 
    50  * 
    51  * 
    52  * SYNOPSIS: 
    53  * 
    54  * double x, y, y0(); 
    55  * 
    56  * y = y0( x ); 
    57  * 
    58  * 
    59  * 
    60  * DESCRIPTION: 
    61  * 
    62  * Returns Bessel function of the second kind, of order 
    63  * zero, of the argument. 
    64  * 
    65  * The domain is divided into the intervals [0, 5] and 
    66  * (5, infinity). In the first interval a rational approximation 
    67  * R(x) is employed to compute 
    68  *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI. 
    69  * Thus a call to j0() is required. 
    70  * 
    71  * In the second interval, the Hankel asymptotic expansion 
    72  * is employed with two rational functions of degree 6/6 
    73  * and 7/7. 
    74  * 
    75  * 
    76  * 
    77  * ACCURACY: 
    78  * 
    79  *  Absolute error, when y0(x) < 1; else relative error: 
    80  * 
    81  * arithmetic   domain     # trials      peak         rms 
    82  *    DEC       0, 30        9400       7.0e-17     7.9e-18 
    83  *    IEEE      0, 30       30000       1.3e-15     1.6e-16 
    84  * 
    85  */ 
    86  
    8746 
    8847/* 
    … …  
    9554 
    9655double j0( double ); 
    97  
    9856double j0(double x) { 
    9957 
    … …  
    291249 
    292250    q = 1.0/x; 
    293     w = sqrtf(q); 
     251    w = sqrt(q); 
    294252 
    295253    p = w * polevl( q, MO, 7); 
    296254    w = q*q; 
    297255    xn = q * polevl( w, PH, 7) - PIO4F; 
    298     p = p * cosf(xn + xx); 
     256    p = p * cos(xn + xx); 
    299257    return(p); 
    300258#endif 

  • sasmodels/models/lib/j1_cephes.c

    rbfef528 re2af2a9  
    3232 *    IEEE      0, 30       30000       2.6e-16     1.1e-16 
    3333 * 
    34  * 
    35  */ 
    36 /*                                                      y1.c 
    37  * 
    38  *      Bessel function of second kind of order one 
    39  * 
    40  * 
    41  * 
    42  * SYNOPSIS: 
    43  * 
    44  * double x, y, y1(); 
    45  * 
    46  * y = y1( x ); 
    47  * 
    48  * 
    49  * 
    50  * DESCRIPTION: 
    51  * 
    52  * Returns Bessel function of the second kind of order one 
    53  * of the argument. 
    54  * 
    55  * The domain is divided into the intervals [0, 8] and 
    56  * (8, infinity). In the first interval a 25 term Chebyshev 
    57  * expansion is used, and a call to j1() is required. 
    58  * In the second, the asymptotic trigonometric representation 
    59  * is employed using two rational functions of degree 5/5. 
    60  * 
    61  * 
    62  * 
    63  * ACCURACY: 
    64  * 
    65  *                      Absolute error: 
    66  * arithmetic   domain      # trials      peak         rms 
    67  *    DEC       0, 30       10000       8.6e-17     1.3e-17 
    68  *    IEEE      0, 30       30000       1.0e-15     1.3e-16 
    69  * 
    70  * (error criterion relative when |y1| > 1). 
    7134 * 
    7235 */ 

  • sasmodels/models/lib/polevl.c

    r3936ad3 re2af2a9  
    5050Direct inquiries to 30 Frost Street, Cambridge, MA 02140 
    5151*/ 
     52 
    5253double polevl( double x, double coef[8], int N ); 
    5354double p1evl( double x, double coef[8], int N ); 

  • sasmodels/models/poly_gauss_coil.py

    r246517d r09b84ed  
    5050""" 
    5151 
    52 from numpy import inf, sqrt, power 
     52from numpy import inf, sqrt, exp, power 
    5353 
    5454name =  "poly_gauss_coil" 
    … …  
    6969def Iq(q, i_zero, radius_gyration, polydispersity): 
    7070    # pylint: disable = missing-docstring 
    71     # need to trap the case of the polydispersity being 1 (ie, monodispersity) 
    7271    u = polydispersity - 1.0 
    73     if polydispersity == 1: 
    74        minusoneonu = -1.0 / u 
    75     else: 
    76        minusoneonu = -1.0 / u 
    7772    z = ((q * radius_gyration) * (q * radius_gyration)) / (1.0 + 2.0 * u) 
    7873    if (q == 0).any(): 
    79        inten = i_zero 
     74        inten = i_zero 
    8075    else: 
    81        inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z)) 
     76    # need to trap the case of the polydispersity being 1 (ie, monodispersity!) 
     77        if polydispersity == 1: 
     78            inten = i_zero * 2.0 * (exp(-z) + z - 1.0 ) / (z * z) 
     79        else: 
     80            minusoneonu = -1.0 / u 
     81            inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z)) 
    8282    return inten 
    83 Iq.vectorized =  True  # Iq accepts an array of q values 
     83#Iq.vectorized =  True  # Iq accepts an array of q values 
    8484 
    8585def Iqxy(qx, qy, *args): 
    … …  
    100100               background = 'background') 
    101101 
     102# these unit test values taken from SasView 3.1.2 
    102103tests =  [ 
    103     [{'scale': 70.0, 'radius_gyration': 75.0, 'polydispersity': 2.0, 'background': 0.0}, 
     104    [{'scale': 1.0, 'i_zero': 70.0, 'radius_gyration': 75.0, 'polydispersity': 2.0, 'background': 0.0}, 
    104105     [0.0106939, 0.469418], [57.6405, 0.169016]], 
    105106    ] 

  • 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 

  • sasmodels/models/adsorbed_layer.py

    rf10bc52 r4f4e0d5  
    66 
    77r""" 
    8 This model describes the scattering from a layer of surfactant or polymer adsorbed on spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for. 
     8This model describes the scattering from a layer of surfactant or polymer adsorbed on large, smooth, notionally spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for. 
    99 
    1010Unlike many other core-shell models, this model does not assume any form for the density distribution of the adsorbed species normal to the interface (cf, a core-shell model normally assumes the density distribution to be a homogeneous step-function). For comparison, if the thickness of a (traditional core-shell like) step function distribution is *t*, the second moment about the mean of the density distribution (ie, the distance of the centre-of-mass of the distribution from the interface), |sigma| =  sqrt((*t* :sup:`2` )/12). 
    … …  
    3434 
    3535description =  """ 
    36     Evaluates the scattering from particles 
     36    Evaluates the scattering from large particles 
    3737    with an adsorbed layer of surfactant or 
    3838    polymer, independent of the form of the 
    … …  
    4242 
    4343#             ["name", "units", default, [lower, upper], "type", "description"], 
    44 parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment"], 
    45               ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount"], 
    46               ["density_poly", "g/cm3", 0.7, [0.0, inf], "", "Polymer density"], 
    47               ["radius", "Ang", 500.0, [0.0, inf], "", "Particle radius"], 
    48               ["volfraction", "None", 0.14, [0.0, inf], "", "Particle vol fraction"], 
    49               ["polymer_sld", "1/Ang^2", 1.5e-06, [-inf, inf], "", "Polymer SLD"], 
    50               ["solvent_sld", "1/Ang^2", 6.3e-06, [-inf, inf], "", "Solvent SLD"]] 
     44parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment of polymer distribution"], 
     45              ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount of polymer"], 
     46              ["density_shell", "g/cm3", 0.7, [0.0, inf], "", "Bulk density of polymer in the shell"], 
     47              ["radius", "Ang", 500.0, [0.0, inf], "", "Core particle radius"], 
     48              ["volfraction", "None", 0.14, [0.0, inf], "", "Core particle volume fraction"], 
     49              ["sld_shell", "1e-6/Ang^2", 1.5, [-inf, inf], "sld", "Polymer shell SLD"], 
     50              ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent SLD"]] 
    5151 
    5252# NB: Scale and Background are implicit parameters on every model 
    53 def Iq(q, second_moment, adsorbed_amount, density_poly, radius,  
    54         volfraction, polymer_sld, solvent_sld): 
     53def Iq(q, second_moment, adsorbed_amount, density_shell, radius,  
     54        volfraction, sld_shell, sld_solvent): 
    5555    # pylint: disable = missing-docstring 
    56     deltarhosqrd =  (polymer_sld - solvent_sld) * (polymer_sld - solvent_sld) 
    57     numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount) 
    58     denominator =  (q * q) * (density_poly * density_poly) * radius 
    59     eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment)) 
    60     #scale by 10^10 for units conversion to cm^-1 
    61     inten =  1.0e+10 * deltarhosqrd * ((numerator / denominator) * eterm) 
     56#    deltarhosqrd =  (sld_shell - sld_solvent) * (sld_shell - sld_solvent) 
     57#    numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount) 
     58#    denominator =  (q * q) * (density_shell * density_shell) * radius 
     59#    eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment)) 
     60#    #scale by 10^-2 for units conversion to cm^-1 
     61#    inten =  1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm) 
     62    aa =  (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell  
     63    bb = q * second_moment 
     64    #scale by 10^-2 for units conversion to cm^-1 
     65    inten =  6.0e-02 * pi * volfraction * aa * aa * exp(-bb * bb) / radius 
    6266    return inten 
    6367Iq.vectorized =  True  # Iq accepts an array of q values 
    … …  
    7175            second_moment = 23.0, 
    7276            adsorbed_amount = 1.9, 
    73             density_poly = 0.7, 
     77            density_shell = 0.7, 
    7478            radius = 500.0, 
    7579            volfraction = 0.14, 
    76             polymer_sld = 1.5e-06, 
    77             solvent_sld = 6.3e-06, 
     80            sld_shell = 1.5, 
     81            sld_solvent = 6.3, 
    7882            background = 0.0) 
    7983 
    … …  
    8286               second_moment = 'second_moment', 
    8387               adsorbed_amount = 'ads_amount', 
    84                density_poly = 'density_poly', 
     88               density_shell = 'density_poly', 
    8589               radius = 'radius_core', 
    8690               volfraction = 'volf_cores', 
    87                polymer_sld = 'sld_poly', 
    88                solvent_sld = 'sld_solv', 
     91               sld_shell = 'sld_poly', 
     92               sld_solvent = 'sld_solv', 
    8993               background = 'background') 
    9094 
    … …  
    9296tests =  [ 
    9397    [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9,  
    94      'density_poly': 0.7, 'radius': 500.0, 'volfraction': 0.14,  
    95      'polymer_sld': 1.5e-06, 'solvent_sld': 6.3e-06, 'background': 0.0}, 
     98     'density_shell': 0.7, 'radius': 500.0, 'volfraction': 0.14,  
     99     'sld_shell': 1.5, 'sld_solvent': 6.3, 'background': 0.0}, 
    96100     [0.0106939, 0.469418], [73.741, 9.65391e-53]], 
    97101    ] 
     102# ADDED by: SMK  ON: 16Mar2016  convert from sasview, check vs SANDRA, 18Mar2016 RKH some edits & renaming 
Note: See TracChangeset for help on using the changeset viewer.