Changes in / [4f4e0d5:eec5fd6] in sasmodels


Ignore:
Files:
1 added
13 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

    r15bd6e7 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 
Note: See TracChangeset for help on using the changeset viewer.