Changeset 17a8c94 in sasmodels for sasmodels


Ignore:
Timestamp:
Mar 31, 2018 9:22:09 PM (6 years ago)
Author:
butler
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
9a7ef88
Parents:
5bc6d21 (diff), b3703f5 (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' into ticket-896

Location:
sasmodels
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r3221de0 rd86f0fc  
    718718    model = core.build_model(model_info, dtype=dtype, platform="ocl") 
    719719    calculator = DirectModel(data, model, cutoff=cutoff) 
    720     engine_type = calculator._model.__class__.__name__.replace('Model','').upper() 
     720    engine_type = calculator._model.__class__.__name__.replace('Model', '').upper() 
    721721    bits = calculator._model.dtype.itemsize*8 
    722722    precision = "fast" if getattr(calculator._model, 'fast', False) else str(bits) 
    … …  
    14911491            vmin, vmax = limits 
    14921492            self.limits = vmax*1e-7, 1.3*vmax 
    1493             import pylab; pylab.clf() 
     1493            import pylab 
     1494            pylab.clf() 
    14941495            plot_models(self.opts, result, limits=self.limits) 
    14951496 

  • sasmodels/data.py

    rf549e37 rd86f0fc  
    706706    else: 
    707707        vmin_scaled, vmax_scaled = vmin, vmax 
    708     nx, ny = len(data.x_bins), len(data.y_bins) 
     708    #nx, ny = len(data.x_bins), len(data.y_bins) 
    709709    x_bins, y_bins, image = _build_matrix(data, plottable) 
    710710    plt.imshow(image, 
    … …  
    772772        if loop >= max_loop:  # this protects never-ending loop 
    773773            break 
    774         image = _fillup_pixels(self, image=image, weights=weights) 
     774        image = _fillup_pixels(image=image, weights=weights) 
    775775        loop += 1 
    776776 
    … …  
    817817    return x_bins, y_bins 
    818818 
    819 def _fillup_pixels(self, image=None, weights=None): 
     819def _fillup_pixels(image=None, weights=None): 
    820820    """ 
    821821    Fill z values of the empty cells of 2d image matrix 

  • sasmodels/generate.py

    r6cbdcd4 rd86f0fc  
    169169 
    170170import sys 
    171 from os.path import abspath, dirname, join as joinpath, exists, getmtime 
     171from os import environ 
     172from os.path import abspath, dirname, join as joinpath, exists, getmtime, sep 
    172173import re 
    173174import string 
    … …  
    289290    import loops. 
    290291    """ 
    291     if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 
    292         import os.path 
     292    if info.source and any(lib.startswith('lib/gauss') for lib in info.source): 
    293293        from .gengauss import gengauss 
    294         path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 
    295         if not os.path.exists(path): 
     294        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n) 
     295        if not exists(path): 
    296296            gengauss(n, path) 
    297297        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 
    298                         else lib for lib in info.source] 
     298                       else lib for lib in info.source] 
    299299 
    300300def format_units(units): 
    … …  
    320320                     for w, h in zip(column_widths, PARTABLE_HEADERS)] 
    321321 
    322     sep = " ".join("="*w for w in column_widths) 
     322    underbar = " ".join("="*w for w in column_widths) 
    323323    lines = [ 
    324         sep, 
     324        underbar, 
    325325        " ".join("%-*s" % (w, h) 
    326326                 for w, h in zip(column_widths, PARTABLE_HEADERS)), 
    327         sep, 
     327        underbar, 
    328328        ] 
    329329    for p in pars: 
    … …  
    334334            "%*g" % (column_widths[3], p.default), 
    335335            ])) 
    336     lines.append(sep) 
     336    lines.append(underbar) 
    337337    return "\n".join(lines) 
    338338 
    … …  
    612612    """ 
    613613    spaces = " "*depth 
    614     sep = "\n" + spaces 
    615     return spaces + sep.join(s.split("\n")) 
     614    interline_separator = "\n" + spaces 
     615    return spaces + interline_separator.join(s.split("\n")) 
    616616 
    617617 
    … …  
    619619def load_template(filename): 
    620620    # type: (str) -> str 
     621    """ 
     622    Load template file from sasmodels resource directory. 
     623    """ 
    621624    path = joinpath(DATA_PATH, filename) 
    622625    mtime = getmtime(path) 
    … …  
    900903        kernel_module = load_custom_kernel_module(model_name) 
    901904    else: 
    902         from sasmodels import models 
    903         __import__('sasmodels.models.'+model_name) 
    904         kernel_module = getattr(models, model_name, None) 
     905        try: 
     906            from sasmodels import models 
     907            __import__('sasmodels.models.'+model_name) 
     908            kernel_module = getattr(models, model_name, None) 
     909        except ImportError: 
     910            # If the model isn't a built in model, try the plugin directory 
     911            plugin_path = environ.get('SAS_MODELPATH', None) 
     912            if plugin_path is not None: 
     913                file_name = model_name.split(sep)[-1] 
     914                model_name = plugin_path + sep + file_name + ".py" 
     915                kernel_module = load_custom_kernel_module(model_name) 
     916            else: 
     917                raise 
    905918    return kernel_module 
    906919 

  • sasmodels/guyou.py

    r0d5a655 rb3703f5  
    3131# 
    3232# 2017-11-01 Paul Kienzle 
    33 # * converted to python, using degrees rather than radians 
     33# * converted to python, with degrees rather than radians 
     34""" 
     35Convert between latitude-longitude and Guyou map coordinates. 
     36""" 
     37 
    3438from __future__ import division, print_function 
    3539 
    3640import numpy as np 
    37 from numpy import sqrt, pi, tan, cos, sin, log, exp, arctan2 as atan2, sign, radians, degrees 
     41from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees 
     42from numpy import sinh, arctan as atan 
     43 
     44# scipy version of special functions 
     45from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 
    3846 
    3947_ = """ 
    … …  
    5967""" 
    6068 
    61 # scipy version of special functions 
    62 from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 
    63 from numpy import sinh, sign, arctan as atan 
    64  
    6569def ellipticJi(u, v, m): 
    6670    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 
    6771    u, v, m = np.broadcast_arrays(u, v, m) 
    68     result = np.empty_like([u,u,u], 'D') 
    69     real = v==0 
    70     imag = u==0 
     72    result = np.empty_like([u, u, u], 'D') 
     73    real = (v == 0) 
     74    imag = (u == 0) 
    7175    mixed = ~(real|imag) 
    7276    result[:, real] = _ellipticJi_real(u[real], m[real]) 
    7377    result[:, imag] = _ellipticJi_imag(v[imag], m[imag]) 
    7478    result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed]) 
    75     return result[0,:] if scalar else result 
     79    return result[0, :] if scalar else result 
    7680 
    7781def _ellipticJi_real(u, m): 
    … …  
    104108        result = np.empty_like(phi, 'D') 
    105109        index = (phi == 0) 
    106         result[index] = ellipticF(atan(sinh(abs(phi[index]))), 1-m[index]) * sign(psi[index]) 
     110        result[index] = ellipticF(atan(sinh(abs(phi[index]))), 
     111                                  1-m[index]) * sign(psi[index]) 
    107112        result[~index] = ellipticFi(phi[~index], psi[~index], m[~index]) 
    108113        return result.reshape(1)[0] if scalar else result 
    … …  
    117122    cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2 
    118123    re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi) 
    119     im = ellipticF(atan(sqrt(np.maximum(0,(cotlambda2 / cotphi2 - 1) / m))), 1 - m) * sign(psi) 
     124    im = ellipticF(atan(sqrt(np.maximum(0, (cotlambda2 / cotphi2 - 1) / m))), 
     125                   1 - m) * sign(psi) 
    120126    return re + 1j*im 
    121127 
    122 sqrt2 = sqrt(2) 
     128SQRT2 = sqrt(2) 
    123129 
    124130# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion 
    … …  
    127133# K = 3.165103454447431823666270142140819753058976299237578486994... 
    128134def guyou(lam, phi): 
     135    """Transform from (latitude, longitude) to point (x, y)""" 
    129136    # [PAK] wrap into [-pi/2, pi/2] radians 
    130137    x, y = np.asarray(lam), np.asarray(phi) 
    … …  
    135142 
    136143    # Compute constant K 
    137     cos_u = (sqrt2 - 1) / (sqrt2 + 1) 
     144    cos_u = (SQRT2 - 1) / (SQRT2 + 1) 
    138145    sinsq_u = 1 - cos_u**2 
    139146    K = ellipticF(pi/2, sinsq_u) 
    … …  
    144151    at = atan(r * (cos(lam) - 1j*sin(lam))) 
    145152    t = ellipticFi(at.real, at.imag, sinsq_u) 
    146     x, y = (-t.imag, sign(phi + (phi==0))*(0.5 * K - t.real)) 
     153    x, y = (-t.imag, sign(phi + (phi == 0))*(0.5 * K - t.real)) 
    147154 
    148155    # [PAK] convert to degrees, and return to original tile 
    … …  
    150157 
    151158def guyou_invert(x, y): 
     159    """Transform from point (x, y) on plot to (latitude, longitude)""" 
    152160    # [PAK] wrap into [-pi/2, pi/2] radians 
    153161    x, y = np.asarray(x), np.asarray(y) 
    … …  
    158166 
    159167    # compute constant K 
    160     cos_u = (sqrt2 - 1) / (sqrt2 + 1) 
     168    cos_u = (SQRT2 - 1) / (SQRT2 + 1) 
    161169    sinsq_u = 1 - cos_u**2 
    162170    K = ellipticF(pi/2, sinsq_u) 
    … …  
    174182 
    175183def plot_grid(): 
     184    """Plot the latitude-longitude grid for Guyou transform""" 
    176185    import matplotlib.pyplot as plt 
    177186    from numpy import linspace 
    … …  
    186195        plt.plot(x, y, 'g') 
    187196 
    188     for long in range(-limit, limit+1, step): 
    189         x, y = guyou(scale*long, scale*long_line) 
     197    for longitude in range(-limit, limit+1, step): 
     198        x, y = guyou(scale*longitude, scale*long_line) 
    190199        plt.plot(x, y, 'b') 
    191200    #plt.xlabel('longitude') 
    192201    plt.ylabel('latitude') 
     202    plt.title('forward transform') 
    193203 
    194204    plt.subplot(212) 
    … …  
    202212    plt.xlabel('longitude') 
    203213    plt.ylabel('latitude') 
     214    plt.title('inverse transform') 
     215 
     216def main(): 
     217    """Show the Guyou transformation""" 
     218    plot_grid() 
     219    import matplotlib.pyplot as plt 
     220    plt.show() 
    204221 
    205222if __name__ == "__main__": 
    206     plot_grid() 
    207     import matplotlib.pyplot as plt; plt.show() 
    208  
     223    main() 
    209224 
    210225_ = """ 

  • sasmodels/jitter.py

    r199bd07 rb3703f5  
    88from __future__ import division, print_function 
    99 
    10 import sys, os 
    11 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
    12  
    1310import argparse 
    1411 
    … …  
    1916 
    2017import matplotlib.pyplot as plt 
    21 from matplotlib.widgets import Slider, CheckButtons 
    22 from matplotlib import cm 
     18from matplotlib.widgets import Slider 
    2319import numpy as np 
    2420from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    2521 
    26 def draw_beam(ax, view=(0, 0)): 
     22def draw_beam(axes, view=(0, 0)): 
    2723    """ 
    2824    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
    2925    """ 
    30     #ax.plot([0,0],[0,0],[1,-1]) 
    31     #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
     26    #axes.plot([0,0],[0,0],[1,-1]) 
     27    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    3228 
    3329    steps = 25 
    … …  
    4642    x, y, z = [v.reshape(shape) for v in points] 
    4743 
    48     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    49  
    50 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     44    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
     45 
     46def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): 
    5147    """Draw an ellipsoid.""" 
    52     a,b,c = size 
     48    a, b, c = size 
    5349    u = np.linspace(0, 2 * np.pi, steps) 
    5450    v = np.linspace(0, np.pi, steps) 
    … …  
    5854    x, y, z = transform_xyz(view, jitter, x, y, z) 
    5955 
    60     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    61  
    62     draw_labels(ax, view, jitter, [ 
    63          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    64          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    65          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    66          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    67          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    68          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     56    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     57 
     58    draw_labels(axes, view, jitter, [ 
     59        ('c+', [+0, +0, +c], [+1, +0, +0]), 
     60        ('c-', [+0, +0, -c], [+0, +0, -1]), 
     61        ('a+', [+a, +0, +0], [+0, +0, +1]), 
     62        ('a-', [-a, +0, +0], [+0, +0, -1]), 
     63        ('b+', [+0, +b, +0], [-1, +0, +0]), 
     64        ('b-', [+0, -b, +0], [-1, +0, +0]), 
    6965    ]) 
    7066 
    71 def draw_sc(ax, size, view, jitter, steps=None, alpha=1): 
     67def draw_sc(axes, size, view, jitter, steps=None, alpha=1): 
     68    """Draw points for simple cubic paracrystal""" 
    7269    atoms = _build_sc() 
    73     _draw_crystal(ax, size, view, jitter, atoms=atoms) 
    74  
    75 def draw_fcc(ax, size, view, jitter, steps=None, alpha=1): 
     70    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     71 
     72def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): 
     73    """Draw points for face-centered cubic paracrystal""" 
    7674    # Build the simple cubic crystal 
    7775    atoms = _build_sc() 
    … …  
    8583    # y and z planes can be generated by substituting x for y and z respectively 
    8684    atoms.extend(zip(x+y+z, y+z+x, z+x+y)) 
    87     _draw_crystal(ax, size, view, jitter, atoms=atoms) 
    88  
    89 def draw_bcc(ax, size, view, jitter, steps=None, alpha=1): 
     85    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     86 
     87def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): 
     88    """Draw points for body-centered cubic paracrystal""" 
    9089    # Build the simple cubic crystal 
    9190    atoms = _build_sc() 
    … …  
    9897    ) 
    9998    atoms.extend(zip(x, y, z)) 
    100     _draw_crystal(ax, size, view, jitter, atoms=atoms) 
    101  
    102 def _draw_crystal(ax, size, view, jitter, steps=None, alpha=1, atoms=None): 
     99    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     100 
     101def _draw_crystal(axes, size, view, jitter, atoms=None): 
    103102    atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd') 
    104103    x, y, z = atoms*size[:, None] 
    105104    x, y, z = transform_xyz(view, jitter, x, y, z) 
    106     ax.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o') 
    107     ax.scatter(x[1:], y[1:], z[1:], c='r', marker='o') 
     105    axes.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o') 
     106    axes.scatter(x[1:], y[1:], z[1:], c='r', marker='o') 
    108107 
    109108def _build_sc(): 
    … …  
    124123    return atoms 
    125124 
    126 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     125def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): 
    127126    """Draw a parallelepiped.""" 
    128127    a, b, c = size 
    129     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    130     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    131     z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     128    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     129    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     130    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    132131    tri = np.array([ 
    133132        # counter clockwise triangles 
    134133        # z: up/down, x: right/left, y: front/back 
    135         [0,1,2], [3,2,1], # top face 
    136         [6,5,4], [5,6,7], # bottom face 
    137         [0,2,6], [6,4,0], # right face 
    138         [1,5,7], [7,3,1], # left face 
    139         [2,3,6], [7,6,3], # front face 
    140         [4,1,0], [5,1,4], # back face 
     134        [0, 1, 2], [3, 2, 1], # top face 
     135        [6, 5, 4], [5, 6, 7], # bottom face 
     136        [0, 2, 6], [6, 4, 0], # right face 
     137        [1, 5, 7], [7, 3, 1], # left face 
     138        [2, 3, 6], [7, 6, 3], # front face 
     139        [4, 1, 0], [5, 1, 4], # back face 
    141140    ]) 
    142141 
    143142    x, y, z = transform_xyz(view, jitter, x, y, z) 
    144     ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
     143    axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    145144 
    146145    # Draw pink face on box. 
    … …  
    149148    # rotate that face. 
    150149    if 1: 
    151         x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    152         y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    153         z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     150        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     151        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     152        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    154153        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
    155         ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha) 
    156  
    157     draw_labels(ax, view, jitter, [ 
    158          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    159          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    160          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    161          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    162          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    163          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     154        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) 
     155 
     156    draw_labels(axes, view, jitter, [ 
     157        ('c+', [+0, +0, +c], [+1, +0, +0]), 
     158        ('c-', [+0, +0, -c], [+0, +0, -1]), 
     159        ('a+', [+a, +0, +0], [+0, +0, +1]), 
     160        ('a-', [-a, +0, +0], [+0, +0, -1]), 
     161        ('b+', [+0, +b, +0], [-1, +0, +0]), 
     162        ('b-', [+0, -b, +0], [-1, +0, +0]), 
    164163    ]) 
    165164 
    166 def draw_sphere(ax, radius=10., steps=100): 
     165def draw_sphere(axes, radius=10., steps=100): 
    167166    """Draw a sphere""" 
    168167    u = np.linspace(0, 2 * np.pi, steps) 
    … …  
    172171    y = radius * np.outer(np.sin(u), np.sin(v)) 
    173172    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    174     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    175  
    176 def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
     173    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
     174 
     175def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
    177176                draw_shape=draw_parallelepiped): 
    178177    """ 
    … …  
    187186    cloud = [ 
    188187        [-1, -1, -1], 
    189         [-1, -1,  0], 
    190         [-1, -1,  1], 
    191         [-1,  0, -1], 
    192         [-1,  0,  0], 
    193         [-1,  0,  1], 
    194         [-1,  1, -1], 
    195         [-1,  1,  0], 
    196         [-1,  1,  1], 
    197         [ 0, -1, -1], 
    198         [ 0, -1,  0], 
    199         [ 0, -1,  1], 
    200         [ 0,  0, -1], 
    201         [ 0,  0,  0], 
    202         [ 0,  0,  1], 
    203         [ 0,  1, -1], 
    204         [ 0,  1,  0], 
    205         [ 0,  1,  1], 
    206         [ 1, -1, -1], 
    207         [ 1, -1,  0], 
    208         [ 1, -1,  1], 
    209         [ 1,  0, -1], 
    210         [ 1,  0,  0], 
    211         [ 1,  0,  1], 
    212         [ 1,  1, -1], 
    213         [ 1,  1,  0], 
    214         [ 1,  1,  1], 
     188        [-1, -1, +0], 
     189        [-1, -1, +1], 
     190        [-1, +0, -1], 
     191        [-1, +0, +0], 
     192        [-1, +0, +1], 
     193        [-1, +1, -1], 
     194        [-1, +1, +0], 
     195        [-1, +1, +1], 
     196        [+0, -1, -1], 
     197        [+0, -1, +0], 
     198        [+0, -1, +1], 
     199        [+0, +0, -1], 
     200        [+0, +0, +0], 
     201        [+0, +0, +1], 
     202        [+0, +1, -1], 
     203        [+0, +1, +0], 
     204        [+0, +1, +1], 
     205        [+1, -1, -1], 
     206        [+1, -1, +0], 
     207        [+1, -1, +1], 
     208        [+1, +0, -1], 
     209        [+1, +0, +0], 
     210        [+1, +0, +1], 
     211        [+1, +1, -1], 
     212        [+1, +1, +0], 
     213        [+1, +1, +1], 
    215214    ] 
    216215    dtheta, dphi, dpsi = jitter 
    … …  
    221220    if dpsi == 0: 
    222221        cloud = [v for v in cloud if v[2] == 0] 
    223     draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     222    draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    224223    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    225224    for point in cloud: 
    226225        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
    227         draw_shape(ax, size, view, delta, alpha=0.8) 
     226        draw_shape(axes, size, view, delta, alpha=0.8) 
    228227    for v in 'xyz': 
    229228        a, b, c = size 
    230         lim = np.sqrt(a**2+b**2+c**2) 
    231         getattr(ax, 'set_'+v+'lim')([-lim, lim]) 
    232         getattr(ax, v+'axis').label.set_text(v) 
     229        lim = np.sqrt(a**2 + b**2 + c**2) 
     230        getattr(axes, 'set_'+v+'lim')([-lim, lim]) 
     231        getattr(axes, v+'axis').label.set_text(v) 
    233232 
    234233PROJECTIONS = [ 
    … …  
    238237    'azimuthal_equal_area', 
    239238] 
    240 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
     239def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
    241240              projection='equirectangular'): 
    242241    """ 
    … …  
    297296        Should allow free movement in theta, but phi is distorted. 
    298297    """ 
    299     theta, phi, psi = view 
    300     dtheta, dphi, dpsi = jitter 
    301  
    302     t = np.linspace(-1, 1, n) 
    303     weights = np.ones_like(t) 
     298    # TODO: try Kent distribution instead of a gaussian warped by projection 
     299 
     300    dist_x = np.linspace(-1, 1, n) 
     301    weights = np.ones_like(dist_x) 
    304302    if dist == 'gaussian': 
    305         t *= 3 
    306         weights = exp(-0.5*t**2) 
     303        dist_x *= 3 
     304        weights = exp(-0.5*dist_x**2) 
    307305    elif dist == 'rectangle': 
    308306        # Note: uses sasmodels ridiculous definition of rectangle width 
    309         t *= sqrt(3) 
     307        dist_x *= sqrt(3) 
    310308    elif dist == 'uniform': 
    311309        pass 
    … …  
    314312 
    315313    if projection == 'equirectangular':  #define PROJECTION 1 
    316         def rotate(theta_i, phi_j): 
     314        def _rotate(theta_i, phi_j): 
    317315            return Rx(phi_j)*Ry(theta_i) 
    318         def weight(theta_i, phi_j, wi, wj): 
    319             return wi*wj*abs(cos(radians(theta_i))) 
     316        def _weight(theta_i, phi_j, w_i, w_j): 
     317            return w_i*w_j*abs(cos(radians(theta_i))) 
    320318    elif projection == 'sinusoidal':  #define PROJECTION 2 
    321         def rotate(theta_i, phi_j): 
     319        def _rotate(theta_i, phi_j): 
    322320            latitude = theta_i 
    323321            scale = cos(radians(latitude)) 
    … …  
    325323            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    326324            return Rx(longitude)*Ry(latitude) 
    327         def weight(theta_i, phi_j, wi, wj): 
     325        def _weight(theta_i, phi_j, w_i, w_j): 
    328326            latitude = theta_i 
    329327            scale = cos(radians(latitude)) 
    330             w = 1 if abs(phi_j) < abs(scale)*180 else 0 
    331             return w*wi*wj 
     328            active = 1 if abs(phi_j) < abs(scale)*180 else 0 
     329            return active*w_i*w_j 
    332330    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    333         def rotate(theta_i, phi_j): 
    334             from guyou import guyou_invert 
     331        def _rotate(theta_i, phi_j): 
     332            from .guyou import guyou_invert 
    335333            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    336334            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    337335            return Rx(longitude[0])*Ry(latitude[0]) 
    338         def weight(theta_i, phi_j, wi, wj): 
    339             return wi*wj 
     336        def _weight(theta_i, phi_j, w_i, w_j): 
     337            return w_i*w_j 
    340338    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    341         def rotate(theta_i, phi_j): 
     339        def _rotate(theta_i, phi_j): 
    342340            latitude = sqrt(theta_i**2 + phi_j**2) 
    343341            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    344342            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    345343            return Rz(longitude)*Ry(latitude) 
    346         def weight(theta_i, phi_j, wi, wj): 
     344        def _weight(theta_i, phi_j, w_i, w_j): 
    347345            # Weighting for each point comes from the integral: 
    348346            #     \int\int I(q, lat, log) sin(lat) dlat dlog 
    … …  
    374372            # the entire sphere, and treats theta and phi identically. 
    375373            latitude = sqrt(theta_i**2 + phi_j**2) 
    376             w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    377             return w*wi*wj if latitude < 180 else 0 
     374            weight = sin(radians(latitude))/latitude if latitude != 0 else 1 
     375            return weight*w_i*w_j if latitude < 180 else 0 
    378376    elif projection == 'azimuthal_equal_area': 
    379         def rotate(theta_i, phi_j): 
    380             R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    381             latitude = 180-degrees(2*np.arccos(R)) 
     377        def _rotate(theta_i, phi_j): 
     378            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
     379            latitude = 180-degrees(2*np.arccos(radius)) 
    382380            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    383381            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    384382            return Rz(longitude)*Ry(latitude) 
    385         def weight(theta_i, phi_j, wi, wj): 
     383        def _weight(theta_i, phi_j, w_i, w_j): 
    386384            latitude = sqrt(theta_i**2 + phi_j**2) 
    387             w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    388             return w*wi*wj if latitude < 180 else 0 
     385            weight = sin(radians(latitude))/latitude if latitude != 0 else 1 
     386            return weight*w_i*w_j if latitude < 180 else 0 
    389387    else: 
    390388        raise ValueError("unknown projection %r"%projection) 
    391389 
    392390    # mesh in theta, phi formed by rotating z 
     391    dtheta, dphi, dpsi = jitter 
    393392    z = np.matrix([[0], [0], [radius]]) 
    394     points = np.hstack([rotate(theta_i, phi_j)*z 
    395                         for theta_i in dtheta*t 
    396                         for phi_j in dphi*t]) 
    397     # select just the active points (i.e., those with phi < 180 
    398     w = np.array([weight(theta_i, phi_j, wi, wj) 
    399                   for wi, theta_i in zip(weights, dtheta*t) 
    400                   for wj, phi_j in zip(weights, dphi*t)]) 
    401     #print(max(w), min(w), min(w[w>0])) 
    402     points = points[:, w>0] 
    403     w = w[w>0] 
    404     w /= max(w) 
    405  
    406     if 0: # Kent distribution 
    407         points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t]) 
    408         xp, yp, zp = [np.array(v).flatten() for v in points] 
    409         kappa = max(1e6, radians(dtheta)/(2*pi)) 
    410         beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa 
    411         w = exp(kappa*zp) #+ beta*(xp**2 + yp**2) 
    412         print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w)) 
    413         #w /= abs(cos(radians( 
    414         #w /= sum(w) 
     393    points = np.hstack([_rotate(theta_i, phi_j)*z 
     394                        for theta_i in dtheta*dist_x 
     395                        for phi_j in dphi*dist_x]) 
     396    dist_w = np.array([_weight(theta_i, phi_j, w_i, w_j) 
     397                       for w_i, theta_i in zip(weights, dtheta*dist_x) 
     398                       for w_j, phi_j in zip(weights, dphi*dist_x)]) 
     399    #print(max(dist_w), min(dist_w), min(dist_w[dist_w > 0])) 
     400    points = points[:, dist_w > 0] 
     401    dist_w = dist_w[dist_w > 0] 
     402    dist_w /= max(dist_w) 
    415403 
    416404    # rotate relative to beam 
    … …  
    419407    x, y, z = [np.array(v).flatten() for v in points] 
    420408    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) 
    421     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    422  
    423 def draw_labels(ax, view, jitter, text): 
     409    axes.scatter(x, y, z, c=dist_w, marker='o', vmin=0., vmax=1.) 
     410 
     411def draw_labels(axes, view, jitter, text): 
    424412    """ 
    425413    Draw text at a particular location. 
    … …  
    435423    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
    436424        zdir = np.asarray(zdir).flatten() 
    437         ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     425        axes.text(p[0], p[1], p[2], label, zdir=zdir) 
    438426 
    439427# Definition of rotation matrices comes from wikipedia: 
    … …  
    441429def Rx(angle): 
    442430    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    443     a = radians(angle) 
    444     R = [[1, 0, 0], 
    445          [0, +cos(a), -sin(a)], 
    446          [0, +sin(a), +cos(a)]] 
    447     return np.matrix(R) 
     431    angle = radians(angle) 
     432    rot = [[1, 0, 0], 
     433           [0, +cos(angle), -sin(angle)], 
     434           [0, +sin(angle), +cos(angle)]] 
     435    return np.matrix(rot) 
    448436 
    449437def Ry(angle): 
    450438    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    451     a = radians(angle) 
    452     R = [[+cos(a), 0, +sin(a)], 
    453          [0, 1, 0], 
    454          [-sin(a), 0, +cos(a)]] 
    455     return np.matrix(R) 
     439    angle = radians(angle) 
     440    rot = [[+cos(angle), 0, +sin(angle)], 
     441           [0, 1, 0], 
     442           [-sin(angle), 0, +cos(angle)]] 
     443    return np.matrix(rot) 
    456444 
    457445def Rz(angle): 
    458446    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    459     a = radians(angle) 
    460     R = [[+cos(a), -sin(a), 0], 
    461          [+sin(a), +cos(a), 0], 
    462          [0, 0, 1]] 
    463     return np.matrix(R) 
     447    angle = radians(angle) 
     448    rot = [[+cos(angle), -sin(angle), 0], 
     449           [+sin(angle), +cos(angle), 0], 
     450           [0, 0, 1]] 
     451    return np.matrix(rot) 
    464452 
    465453def transform_xyz(view, jitter, x, y, z): 
    … …  
    469457    x, y, z = [np.asarray(v) for v in (x, y, z)] 
    470458    shape = x.shape 
    471     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     459    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
    472460    points = apply_jitter(jitter, points) 
    473461    points = orient_relative_to_beam(view, points) 
    … …  
    526514        return data[offset], data[-1] 
    527515 
    528 def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): 
     516def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): 
    529517    """ 
    530518    Plot the scattering for the particular view. 
    531519 
    532     *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes 
     520    *calculator* is returned from :func:`build_model`.  *axes* are the 3D axes 
    533521    on which the data will be plotted.  *view* and *jitter* are the current 
    534522    orientation and orientation dispersity.  *dist* is one of the sasmodels 
    … …  
    543531    theta, phi, psi = view 
    544532    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
    545     theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     533    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd > 0, phi_pd > 0, psi_pd > 0)] 
    546534    ## increase pd_n for testing jitter integration rather than simple viz 
    547535    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
    … …  
    571559    if 0: 
    572560        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    573         level[level<0] = 0 
     561        level[level < 0] = 0 
    574562        colors = plt.get_cmap()(level) 
    575         ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     563        axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
    576564    elif 1: 
    577         ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
    578                     levels=np.linspace(vmin, vmax, 24)) 
     565        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     566                      levels=np.linspace(vmin, vmax, 24)) 
    579567    else: 
    580         ax.pcolormesh(qx, qy, Iqxy) 
     568        axes.pcolormesh(qx, qy, Iqxy) 
    581569 
    582570def build_model(model_name, n=150, qmax=0.5, **pars): 
    … …  
    595583    for details. 
    596584    """ 
    597     from sasmodels.core import load_model_info, build_model 
     585    from sasmodels.core import load_model_info, build_model as build_sasmodel 
    598586    from sasmodels.data import empty_data2D 
    599587    from sasmodels.direct_model import DirectModel 
    600588 
    601589    model_info = load_model_info(model_name) 
    602     model = build_model(model_info) #, dtype='double!') 
     590    model = build_sasmodel(model_info) #, dtype='double!') 
    603591    q = np.linspace(-qmax, qmax, n) 
    604592    data = empty_data2D(q, q) 
    … …  
    618606    return calculator 
    619607 
    620 def select_calculator(model_name, n=150, size=(10,40,100)): 
     608def select_calculator(model_name, n=150, size=(10, 40, 100)): 
    621609    """ 
    622610    Create a model calculator for the given shape. 
    … …  
    641629        radius = 0.5*c 
    642630        calculator = build_model('sc_paracrystal', n=n, dnn=dnn, 
    643                                   d_factor=d_factor, radius=(1-d_factor)*radius, 
    644                                   background=0) 
     631                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     632                                 background=0) 
    645633    elif model_name == 'fcc_paracrystal': 
    646634        a = b = c 
    … …  
    650638        radius = sqrt(2)/4 * c 
    651639        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn, 
    652                                   d_factor=d_factor, radius=(1-d_factor)*radius, 
    653                                   background=0) 
     640                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     641                                 background=0) 
    654642    elif model_name == 'bcc_paracrystal': 
    655643        a = b = c 
    … …  
    659647        radius = sqrt(3)/2 * c 
    660648        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn, 
    661                                   d_factor=d_factor, radius=(1-d_factor)*radius, 
    662                                   background=0) 
     649                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     650                                 background=0) 
    663651    elif model_name == 'cylinder': 
    664652        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
    … …  
    685673    'cylinder', 
    686674    'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
    687  ] 
     675] 
    688676 
    689677DRAW_SHAPES = { 
    … …  
    751739    plt.gcf().canvas.set_window_title(projection) 
    752740    #gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) 
    753     #ax = plt.subplot(gs[0], projection='3d') 
    754     ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
     741    #axes = plt.subplot(gs[0], projection='3d') 
     742    axes = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
    755743    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection 
    756         ax.axis('square') 
     744        axes.axis('square') 
    757745    except Exception: 
    758746        pass 
    … …  
    761749 
    762750    ## add control widgets to plot 
    763     axtheta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    764     axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    765     axpsi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
    766     stheta = Slider(axtheta, 'Theta', -90, 90, valinit=theta) 
    767     sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    768     spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
    769  
    770     axdtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    771     axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    772     axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     751    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
     752    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
     753    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
     754    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
     755    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
     756    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
     757 
     758    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
     759    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
     760    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    773761    # Note: using ridiculous definition of rectangle distribution, whose width 
    774762    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    775763    # the maximum width to 90. 
    776764    dlimit = DIST_LIMITS[dist] 
    777     sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
    778     sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
    779     sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     765    sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
     766    sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     767    sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
    780768 
    781769 
    … …  
    788776        limit = [0, 0, 0.5, 5][dims] 
    789777        jitter = [0 if v < limit else v for v in jitter] 
    790         ax.cla() 
    791         draw_beam(ax, (0, 0)) 
    792         draw_jitter(ax, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
    793         #draw_jitter(ax, view, (0,0,0)) 
    794         draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
    795         draw_scattering(calculator, ax, view, jitter, dist=dist) 
     778        axes.cla() 
     779        draw_beam(axes, (0, 0)) 
     780        draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
     781        #draw_jitter(axes, view, (0,0,0)) 
     782        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) 
     783        draw_scattering(calculator, axes, view, jitter, dist=dist) 
    796784        plt.gcf().canvas.draw() 
    797785 
    798786    ## bind control widgets to view updater 
    799     stheta.on_changed(lambda v: update(v,'theta')) 
     787    stheta.on_changed(lambda v: update(v, 'theta')) 
    800788    sphi.on_changed(lambda v: update(v, 'phi')) 
    801789    spsi.on_changed(lambda v: update(v, 'psi')) 
    … …  
    815803        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    816804        ) 
    817     parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
    818     parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
    819     parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
    820     parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
    821     parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     805    parser.add_argument('-p', '--projection', choices=PROJECTIONS, 
     806                        default=PROJECTIONS[0], 
     807                        help='coordinate projection') 
     808    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
     809                        help='a,b,c lengths') 
     810    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
     811                        default=DISTRIBUTIONS[0], 
     812                        help='jitter distribution') 
     813    parser.add_argument('-m', '--mesh', type=int, default=30, 
     814                        help='#points in theta-phi mesh') 
     815    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 
     816                        help='oriented shape') 
    822817    opts = parser.parse_args() 
    823818    size = tuple(int(v) for v in opts.size.split(',')) 

  • sasmodels/kernel_iq.c

    raadec17 rd86f0fc  
    173173// Apply the rotation matrix returned from qac_rotation to the point (qx,qy), 
    174174// returning R*[qx,qy]' = [qa,qc]' 
    175 static double 
     175static void 
    176176qac_apply( 
    177177    QACRotation *rotation, 
    … …  
    246246// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy), 
    247247// returning R*[qx,qy]' = [qa,qb,qc]' 
    248 static double 
     248static void 
    249249qabc_apply( 
    250250    QABCRotation *rotation, 

  • sasmodels/kernelcl.py

    r6cbdcd4 rd86f0fc  
    151151        if not HAVE_OPENCL: 
    152152            raise RuntimeError("OpenCL startup failed with ***" 
    153                             + OPENCL_ERROR + "***; using C compiler instead") 
     153                               + OPENCL_ERROR + "***; using C compiler instead") 
    154154        reset_environment() 
    155155        if ENV is None: 

  • sasmodels/models/core_shell_cylinder.c

    r108e70e rd86f0fc  
    4848 
    4949 
    50 double Iqac(double qab, double qc, 
     50static double 
     51Iqac(double qab, double qc, 
    5152    double core_sld, 
    5253    double shell_sld, 

  • sasmodels/models/hollow_rectangular_prism.c

    r108e70e rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio, double thickness); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    63{ 
    74    double length_b = length_a * b2a_ratio; 
    … …  
    1613} 
    1714 
    18 double Iq(double q, 
     15static double 
     16Iq(double q, 
    1917    double sld, 
    2018    double solvent_sld, 
    … …  
    8583} 
    8684 
    87 double Iqabc(double qa, double qb, double qc, 
     85static double 
     86Iqabc(double qa, double qb, double qc, 
    8887    double sld, 
    8988    double solvent_sld, 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r74768cb rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    63{ 
    74    double length_b = length_a * b2a_ratio; 
    … …  
    118} 
    129 
    13 double Iq(double q, 
     10static double 
     11Iq(double q, 
    1412    double sld, 
    1513    double solvent_sld, 

  • sasmodels/models/rectangular_prism.c

    r108e70e rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    63{ 
    74    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 
    85} 
    96 
    10 double Iq(double q, 
     7static double 
     8Iq(double q, 
    119    double sld, 
    1210    double solvent_sld, 
    … …  
    7169 
    7270 
    73 double Iqabc(double qa, double qb, double qc, 
     71static double 
     72Iqabc(double qa, double qb, double qc, 
    7473    double sld, 
    7574    double solvent_sld, 

  • sasmodels/multiscat.py

    r49d1f8b8 rb3703f5  
    7373import argparse 
    7474import time 
    75 import os.path 
    7675 
    7776import numpy as np 
    … …  
    8180from sasmodels import core 
    8281from sasmodels import compare 
    83 from sasmodels import resolution2d 
    8482from sasmodels.resolution import Resolution, bin_edges 
    85 from sasmodels.data import empty_data1D, empty_data2D, plot_data 
    8683from sasmodels.direct_model import call_kernel 
    8784import sasmodels.kernelcl 
    … …  
    106103USE_FAST = True  # OpenCL faster, less accurate math 
    107104 
    108 class NumpyCalculator: 
     105class ICalculator: 
     106    """ 
     107    Multiple scattering calculator 
     108    """ 
     109    def fft(self, Iq): 
     110        """ 
     111        Compute the forward FFT for an image, real -> complex. 
     112        """ 
     113        raise NotImplementedError() 
     114 
     115    def ifft(self, Iq): 
     116        """ 
     117        Compute the inverse FFT for an image, complex -> complex. 
     118        """ 
     119        raise NotImplementedError() 
     120 
     121    def mulitple_scattering(self, Iq): 
     122        r""" 
     123        Compute multiple scattering for I(q) given scattering probability p. 
     124 
     125        Given a probability p of scattering with the thickness, the expected 
     126        number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a 
     127        Poisson weighted sum of single, double, triple, etc. scattering patterns. 
     128        The number of patterns used is based on coverage (default 99%). 
     129        """ 
     130        raise NotImplementedError() 
     131 
     132class NumpyCalculator(ICalculator): 
     133    """ 
     134    Multiple scattering calculator using numpy fft. 
     135    """ 
    109136    def __init__(self, dims=None, dtype=PRECISION): 
    110137        self.dtype = dtype 
    111138        self.complex_dtype = np.dtype('F') if dtype == np.dtype('f') else np.dtype('D') 
    112         pass 
    113139 
    114140    def fft(self, Iq): 
    … …  
    127153 
    128154    def multiple_scattering(self, Iq, p, coverage=0.99): 
    129         r""" 
    130         Compute multiple scattering for I(q) given scattering probability p. 
    131  
    132         Given a probability p of scattering with the thickness, the expected 
    133         number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a 
    134         Poisson weighted sum of single, double, triple, etc. scattering patterns. 
    135         The number of patterns used is based on coverage (default 99%). 
    136         """ 
    137155        #t0 = time.time() 
    138156        coeffs = scattering_coeffs(p, coverage) 
    … …  
    140158        scale = np.sum(Iq) 
    141159        frame = _forward_shift(Iq/scale, dtype=self.dtype) 
    142         F = np.fft.fft2(frame) 
    143         F_convolved = F * np.polyval(poly, F) 
    144         frame = np.fft.ifft2(F_convolved) 
     160        fourier_frame = np.fft.fft2(frame) 
     161        convolved = fourier_frame * np.polyval(poly, fourier_frame) 
     162        frame = np.fft.ifft2(convolved) 
    145163        result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
    146164        #print("numpy multiscat time", time.time()-t0) 
    … …  
    173191""" 
    174192 
    175 class OpenclCalculator(NumpyCalculator): 
     193class OpenclCalculator(ICalculator): 
     194    """ 
     195    Multiple scattering calculator using OpenCL via pyfft. 
     196    """ 
    176197    polyval1f = None 
    177198    polyval1d = None 
    … …  
    180201        context = env.get_context(dtype) 
    181202        if dtype == np.dtype('f'): 
    182             if self.polyval1f is None: 
     203            if OpenclCalculator.polyval1f is None: 
    183204                program = sasmodels.kernelcl.compile_model( 
    184205                    context, POLYVAL1_KERNEL, dtype, fast=USE_FAST) 
    … …  
    187208            self.dtype = dtype 
    188209            self.complex_dtype = np.dtype('F') 
    189             self.polyval1 = self.polyval1f 
     210            self.polyval1 = OpenclCalculator.polyval1f 
    190211        else: 
    191             if self.polyval1d is None: 
     212            if OpenclCalculator.polyval1d is None: 
    192213                program = sasmodels.kernelcl.compile_model( 
    193214                    context, POLYVAL1_KERNEL, dtype, fast=False) 
    … …  
    196217            self.dtype = dtype 
    197218            self.complex_type = np.dtype('D') 
    198             self.polyval1 = self.polyval1d 
     219            self.polyval1 = OpenclCalculator.polyval1d 
    199220        self.queue = env.get_queue(dtype) 
    200221        self.plan = pyfft.cl.Plan(dims, queue=self.queue) 
    … …  
    229250        gpu_poly = cl_array.to_device(self.queue, poly) 
    230251        self.plan.execute(gpu_data.data) 
    231         degree, n = poly.shape[0], frame.shape[0]*frame.shape[1] 
     252        degree, data_size= poly.shape[0], frame.shape[0]*frame.shape[1] 
    232253        self.polyval1( 
    233             self.queue, [n], None, 
    234             np.int32(degree), gpu_poly.data, np.int32(n), gpu_data.data) 
     254            self.queue, [data_size], None, 
     255            np.int32(degree), gpu_poly.data, np.int32(data_size), gpu_data.data) 
    235256        self.plan.execute(gpu_data.data, inverse=True) 
    236257        frame = gpu_data.get() 
    … …  
    251272    """ 
    252273    if transform is None: 
    253         nx, ny = Iq.shape 
    254         transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) 
     274        n_x, n_y = Iq.shape 
     275        transform = Calculator(dims=(n_x*2, n_y*2), dtype=dtype) 
    255276    scale = np.sum(Iq) 
    256277    frame = _forward_shift(Iq/scale, dtype=dtype) 
    … …  
    261282 
    262283def scattering_coeffs(p, coverage=0.99): 
     284    r""" 
     285    Return the coefficients of the scattering powers for transmission 
     286    probability *p*.  This is just the corresponding values for the 
     287    Poisson distribution for $\lambda = -\ln(1-p)$ such that 
     288    $\sum_{k = 0 \ldots n} P(k; \lambda)$ is larger than *coverage*. 
     289    """ 
    263290    L = -np.log(1-p) 
    264291    num_scatter = truncated_poisson_invcdf(coverage, L) 
    … …  
    266293    return coeffs 
    267294 
    268 def truncated_poisson_invcdf(p, L): 
     295def truncated_poisson_invcdf(coverage, L): 
    269296    r""" 
    270     Return smallest k such that cdf(k; L) > p from the truncated Poisson 
     297    Return smallest k such that cdf(k; L) > coverage from the truncated Poisson 
    271298    probability excluding k=0 
    272299    """ 
    … …  
    275302    pmf = -np.exp(-L) / np.expm1(-L) 
    276303    k = 0 
    277     while cdf < p: 
     304    while cdf < coverage: 
    278305        k += 1 
    279306        pmf *= L/k 
    … …  
    305332 
    306333class MultipleScattering(Resolution): 
     334    r""" 
     335    Compute multiple scattering using Fourier convolution. 
     336 
     337    The fourier steps are determined by *qmax*, the maximum $q$ value 
     338    desired, *nq* the number of $q$ steps and *window*, the amount 
     339    of padding around the circular convolution.  The $q$ spacing 
     340    will be $\Delta q = 2 q_\mathrm{max} w / n_q$.  If *nq* is not 
     341    given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$. 
     342 
     343    *probability* is related to the expected number of scattering 
     344    events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
     345    hack to allow probability to be a fitted parameter, the "value" 
     346    can be a function that takes no parameters and returns the current 
     347    value of the probability.  *coverage* determines how many scattering 
     348    steps to consider.  The default is 0.99, which sets $n$ such that 
     349    $1 \ldots n$ covers 99% of the Poisson probability mass function. 
     350 
     351    *is2d* is True then 2D scattering is used, otherwise it accepts 
     352    and returns 1D scattering. 
     353 
     354    *resolution* is the resolution function to apply after multiple 
     355    scattering.  If present, then the resolution $q$ vectors will provide 
     356    default values for *qmin*, *qmax* and *nq*. 
     357    """ 
    307358    def __init__(self, qmin=None, qmax=None, nq=None, window=2, 
    308359                 probability=None, coverage=0.99, 
    309360                 is2d=False, resolution=None, 
    310361                 dtype=PRECISION): 
    311         r""" 
    312         Compute multiple scattering using Fourier convolution. 
    313  
    314         The fourier steps are determined by *qmax*, the maximum $q$ value 
    315         desired, *nq* the number of $q$ steps and *window*, the amount 
    316         of padding around the circular convolution.  The $q$ spacing 
    317         will be $\Delta q = 2 q_\mathrm{max} w / n_q$.  If *nq* is not 
    318         given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$. 
    319  
    320         *probability* is related to the expected number of scattering 
    321         events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
    322         hack to allow probability to be a fitted parameter, the "value" 
    323         can be a function that takes no parameters and returns the current 
    324         value of the probability.  *coverage* determines how many scattering 
    325         steps to consider.  The default is 0.99, which sets $n$ such that 
    326         $1 \ldots n$ covers 99% of the Poisson probability mass function. 
    327  
    328         *is2d* is True then 2D scattering is used, otherwise it accepts 
    329         and returns 1D scattering. 
    330  
    331         *resolution* is the resolution function to apply after multiple 
    332         scattering.  If present, then the resolution $q$ vectors will provide 
    333         default values for *qmin*, *qmax* and *nq*. 
    334         """ 
    335362        # Infer qmin, qmax from instrument resolution calculator, if present 
    336363        if resolution is not None: 
    … …  
    424451        # Prepare the multiple scattering calculator (either numpy or OpenCL) 
    425452        self.transform = Calculator((2*nq, 2*nq), dtype=dtype) 
     453 
     454        # Iq and Iqxy will be set during apply 
     455        self.Iq = None # type: np.ndarray 
     456        self.Iqxy = None # type: np.ndarray 
    426457 
    427458    def apply(self, theory): 
    … …  
    471502 
    472503    def radial_profile(self, Iqxy): 
     504        """ 
     505        Compute that radial profile for the given Iqxy grid.  The grid should 
     506        be defined as for 
     507        """ 
    473508        # circular average, no anti-aliasing 
    474509        Iq = np.histogram(self._radius, bins=self._edges, weights=Iqxy)[0]/self._norm 
    … …  
    478513def annular_average(qxy, Iqxy, qbins): 
    479514    """ 
    480     Compute annular average of points at 
     515    Compute annular average of points in *Iqxy* at *qbins*.  The $q_x$, $q_y$ 
     516    coordinates for *Iqxy* are given in *qxy*. 
    481517    """ 
    482518    qxy, Iqxy = qxy.flatten(), Iqxy.flatten() 
    … …  
    513549def parse_pars(model, opts): 
    514550    # type: (ModelInfo, argparse.Namespace) -> Dict[str, float] 
     551    """ 
     552    Parse par=val arguments from the command line. 
     553    """ 
    515554 
    516555    seed = np.random.randint(1000000) if opts.random and opts.seed < 0 else opts.seed 
    … …  
    526565        'is2d': opts.is2d, 
    527566    } 
    528     pars, pars2 = compare.parse_pars(compare_opts) 
     567    # Note: sascomp allows comparison on a pair of models, so ignore the second. 
     568    pars, _ = compare.parse_pars(compare_opts) 
    529569    return pars 
    530570 
    … …  
    535575        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    536576        ) 
    537     parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") 
    538     parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') 
    539     parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q') 
    540     parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window') 
    541     parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample') 
    542     parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed') 
    543     parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed') 
    544     parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') 
    545     parser.add_argument('model', type=str, help='sas model name such as cylinder') 
    546     parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') 
     577    parser.add_argument('-p', '--probability', type=float, default=0.1, 
     578                        help="scattering probability") 
     579    parser.add_argument('-n', '--nq', type=int, default=1024, 
     580                        help='number of mesh points') 
     581    parser.add_argument('-q', '--qmax', type=float, default=0.5, 
     582                        help='max q') 
     583    parser.add_argument('-w', '--window', type=float, default=2.0, 
     584                        help='q calc = q max * window') 
     585    parser.add_argument('-2', '--2d', dest='is2d', action='store_true', 
     586                        help='oriented sample') 
     587    parser.add_argument('-s', '--seed', default=-1, 
     588                        help='random pars with given seed') 
     589    parser.add_argument('-r', '--random', action='store_true', 
     590                        help='random pars with random seed') 
     591    parser.add_argument('-o', '--outfile', type=str, default="", 
     592                        help='random pars with random seed') 
     593    parser.add_argument('model', type=str, 
     594                        help='sas model name such as cylinder') 
     595    parser.add_argument('pars', type=str, nargs='*', 
     596                        help='model parameters such as radius=30') 
    547597    opts = parser.parse_args() 
    548598    assert opts.nq%2 == 0, "require even # points" 
    … …  
    592642            plotxy((res._q_steps, res._q_steps), res.Iqxy+background) 
    593643            pylab.title("total scattering for p=%g" % probability) 
     644            if res.resolution is not None: 
     645                pylab.figure() 
     646                plotxy((res._q_steps, res._q_steps), result) 
     647                pylab.title("total scattering with resolution") 
    594648    else: 
    595649        q = res._q 
    … …  
    609663            # Plot 1D pattern for partial scattering 
    610664            pylab.loglog(q, res.Iq+background, label="total for p=%g"%probability) 
     665            if res.resolution is not None: 
     666                pylab.loglog(q, result, label="total with dQ") 
    611667            #new_annulus = annular_average(res._radius, res.Iqxy, res._edges) 
    612668            #pylab.loglog(q, new_annulus+background, label="new total for p=%g"%probability) 

  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rdbf1a60  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
     61    // substitute d_cos_alpha for sin_alpha d_alpha 
    6162    double outer_sum = 0; //initialize integral 
    6263    for( int i=0; i<GAUSS_N; i++) { 
    6364        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6465        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    65  
    66         // inner integral (with gauss points), integration limits = 0, pi/2 
    6766        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6867        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     68 
     69        // inner integral (with gauss points), integration limits = 0, 1 
     70        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    6971        double inner_sum = 0.0; 
    7072        for(int j=0; j<GAUSS_N; j++) { 
    71             const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     73            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    7274            double sin_beta, cos_beta; 
    73             SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     75            SINCOS(M_PI_2*u, sin_beta, cos_beta); 
    7476            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7577            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
    … …  
    9193            inner_sum += GAUSS_W[j] * f * f; 
    9294        } 
     95        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9396        inner_sum *= 0.5; 
    9497        // now sum up the outer integral 
    9598        outer_sum += GAUSS_W[i] * inner_sum; 
    9699    } 
     100    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    97101    outer_sum *= 0.5; 
    98102 

  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 r5bc6d21  
    1111.. math:: 
    1212 
    13     I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
     13    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle  
     14    + \text{background} 
    1415 
    1516where $\langle \ldots \rangle$ is an average over all possible orientations 
    16 of the rectangular solid. 
    17  
    18 The function calculated is the form factor of the rectangular solid below. 
     17of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
     18pulled out of the form factor term due to the multiple slds in the model. 
     19 
    1920The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
    2021$A < B < C$. 
    2122 
    22 .. image:: img/core_shell_parallelepiped_geometry.jpg 
     23.. figure:: img/parallelepiped_geometry.jpg 
     24 
     25   Core of the core shell parallelepiped with the corresponding definition 
     26   of sides. 
     27 
    2328 
    2429There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    2530(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
    26 $(=t_C)$ faces. The projection in the $AB$ plane is then 
    27  
    28 .. image:: img/core_shell_parallelepiped_projection.jpg 
    29  
    30 The volume of the solid is 
     31$(=t_C)$ faces. The projection in the $AB$ plane is 
     32 
     33.. figure:: img/core_shell_parallelepiped_projection.jpg 
     34 
     35   AB cut through the core-shell parallelipiped showing the cross secion of 
     36   four of the six shell slabs. As can be seen, this model leaves **"gaps"** 
     37   at the corners of the solid. 
     38 
     39 
     40The total volume of the solid is thus given as 
    3141 
    3242.. math:: 
    3343 
    3444    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    35  
    36 **meaning that there are "gaps" at the corners of the solid.** 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-shell intensity being calculated as the square of the sum of the 
    40 amplitudes of the core and the slabs on the edges. 
    41  
    42 the scattering amplitude is computed for a particular orientation of the 
    43 core-shell parallelepiped with respect to the scattering vector and then 
    44 averaged over all possible orientations, where $\alpha$ is the angle between 
    45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    46 the angle between projection of the particle in the $xy$ detector plane 
    47 and the $y$ axis. 
    48  
    49 .. math:: 
    50  
    51     F(Q) 
     48amplitudes of the core and the slabs on the edges. The scattering amplitude is 
     49computed for a particular orientation of the core-shell parallelepiped with 
     50respect to the scattering vector and then averaged over all possible 
     51orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
     52of the parallelepiped, and $\beta$ is the angle between the projection of the 
     53particle in the $xy$ detector plane and the $y$ axis. 
     54 
     55.. math:: 
     56 
     57    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha 
     58    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} 
     59 
     60and 
     61 
     62.. math:: 
     63 
     64    F(q,\alpha,\beta) 
    5265    &= (\rho_\text{core}-\rho_\text{solvent}) 
    5366       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    5467    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    55         \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     68        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
    5669    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    5770        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
    … …  
    6376.. math:: 
    6477 
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     78    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 
    6679 
    6780and 
    … …  
    6982.. math:: 
    7083 
    71     Q_A &= \sin\alpha \sin\beta \\ 
    72     Q_B &= \sin\alpha \cos\beta \\ 
    73     Q_C &= \cos\alpha 
     84    Q_A &= q \sin\alpha \sin\beta \\ 
     85    Q_B &= q \sin\alpha \cos\beta \\ 
     86    Q_C &= q \cos\alpha 
    7487 
    7588 
    7689where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    77 are the scattering length of the parallelepiped core, and the rectangular 
     90are the scattering lengths of the parallelepiped core, and the rectangular 
    7891slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    7992is the scattering length of the solvent. 
     93 
     94.. note::  
     95 
     96   the code actually implements two substitutions: $d(cos\alpha)$ is 
     97   substituted for -$sin\alpha \ d\alpha$ (note that in the 
     98   :ref:`parallelepiped` code this is explicitly implemented with 
     99   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that 
     100   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0 
     101   to $\pi/2$. 
    80102 
    81103FITTING NOTES 
    … …  
    94116based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    95117and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96 to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     118to give an oblate or prolate particle, to give an effective radius 
     119for $S(q)$ when $P(q) * S(q)$ is applied. 
    98120 
    99121For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     122angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further 
    101123details of the calculation and angular dispersions see :ref:`orientation`. 
    102124The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    103125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104126 
    105 For 2d, constraints must be applied during fitting to ensure that the 
    106 inequality $A < B < C$ is not violated, and hence the correct definition 
    107 of angles is preserved. The calculation will not report an error, 
    108 but the results may be not correct. 
     127.. note:: For 2d, constraints must be applied during fitting to ensure that the 
     128   inequality $A < B < C$ is not violated, and hence the correct definition 
     129   of angles is preserved. The calculation will not report an error, 
     130   but the results may be not correct. 
    109131 
    110132.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    113135    Note that rotation $\theta$, initially in the $xz$ plane, is carried 
    114136    out first, then rotation $\phi$ about the $z$ axis, finally rotation 
    115     $\Psi$ is now around the axis of the cylinder. The neutron or X-ray 
     137    $\Psi$ is now around the axis of the particle. The neutron or X-ray 
    116138    beam is along the $z$ axis. 
    117139 
    … …  
    120142    Examples of the angles for oriented core-shell parallelepipeds against the 
    121143    detector plane. 
     144 
     145 
     146Validation 
     147---------- 
     148 
     149Cross-checked against hollow rectangular prism and rectangular prism for equal 
     150thickness overlapping sides, and by Monte Carlo sampling of points within the 
     151shape for non-uniform, non-overlapping sides. 
     152 
    122153 
    123154References 
    … …  
    135166 
    136167* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
     168* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 
    138169* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
    139 * Cross-checked against hollow rectangular prism and rectangular prism for 
    140   equal thickness overlapping sides, and by Monte Carlo sampling of points 
    141   within the shape for non-uniform, non-overlapping sides. 
    142170""" 
    143171 

  • sasmodels/models/parallelepiped.c

    r108e70e rdbf1a60  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
     40        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4041        inner_total *= 0.5; 
    4142 
    … …  
    4344        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4445    } 
     46    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4547    outer_total *= 0.5; 
    4648 

  • sasmodels/models/parallelepiped.py

    ref07e95 r5bc6d21  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 The form factor is normalized by the particle volume. 
    5 For information about polarised and magnetic scattering, see 
    6 the :ref:`magnetism` documentation. 
    7  
    84Definition 
    95---------- 
    106 
    117 This model calculates the scattering from a rectangular parallelepiped 
    12  (\:numref:`parallelepiped-image`\). 
    13  If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
     8 (:numref:`parallelepiped-image`). 
     9 If you need to apply polydispersity, see also :ref:`rectangular-prism`. For 
     10 information about polarised and magnetic scattering, see 
     11the :ref:`magnetism` documentation. 
    1412 
    1513.. _parallelepiped-image: 
    … …  
    2624error, or fixing of some dimensions at expected values, may help. 
    2725 
    28 The 1D scattering intensity $I(q)$ is calculated as: 
     26The form factor is normalized by the particle volume and the 1D scattering 
     27intensity $I(q)$ is then calculated as: 
    2928 
    3029.. Comment by Miguel Gonzalez: 
    … …  
    3938 
    4039    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    41            \left< P(q, \alpha) \right> + \text{background} 
     40           \left< P(q, \alpha, \beta) \right> + \text{background} 
    4241 
    4342where the volume $V = A B C$, the contrast is defined as 
    44 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 
    45 $P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
    46 at an angle $\alpha$ (angle between the long axis C and $\vec q$), 
    47 and the averaging $\left<\ldots\right>$ is applied over all orientations. 
     43$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
     44is the form factor corresponding to a parallelepiped oriented 
     45at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 
     46( the angle between the projection of the particle in the $xy$ detector plane 
     47and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
     48orientations. 
    4849 
    4950Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    50 form factor is given by (Mittelbach and Porod, 1961) 
     51form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
    5152 
    5253.. math:: 
    … …  
    6667    \mu &= qB 
    6768 
    68 The scattering intensity per unit volume is returned in units of |cm^-1|. 
     69where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
     70applied. 
    6971 
    7072NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    … …  
    120122.. math:: 
    121123 
    122     P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 
    123                   \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 
    124                   \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 
     124    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 
     125                   {2}qA\cos\alpha)}\right]^2 
     126                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 
     127                   {2}qB\cos\beta)}\right]^2 
     128                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 
     129                   {2}qC\cos\gamma)}\right]^2 
    125130 
    126131with 
    … …  
    160165---------- 
    161166 
    162 P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
    163  
    164 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     167.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 
     168   14 (1961) 185-211 
     169.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    165170 
    166171Authorship and Verification 
Note: See TracChangeset for help on using the changeset viewer.