Changes in / [352494a:956960a] in sasmodels


Ignore:
Files:
14 added
15 deleted
50 edited

Legend:

Unmodified
Added
Removed

  • .travis.yml

    rbec0e75 r24cd982  
    1 # Test Travis CL 
     1language: python 
    22 
    3 language: python 
    4 python: 
    5   - "2.7" 
     3sudo:  false 
     4 
     5matrix: 
     6  include: 
     7    - os: linux 
     8      env: 
     9        - PY=2.7 
     10        - NUMPYSPEC=numpy 
     11    - os: linux 
     12      env: 
     13        - PY=3 
     14        - NUMPYSPEC=numpy 
     15    - os: osx 
     16      language: generic 
     17      env: 
     18        - PY=2.7 
     19        - NUMPYSPEC=numpy 
     20    - os: osx 
     21      language: generic 
     22      env: 
     23        - PY=3 
     24        - NUMPYSPEC=numpy 
     25 
    626# whitelist 
    727branches: 
    828  only: 
    929    - master 
    10 # command to install dependencies 
    11 virtualenv: 
    12   system_site_packages: true 
     30 
     31addons: 
     32  apt: 
     33    packages: 
     34      opencl-headers 
     35 
    1336before_install: 
    14   - sudo apt-get update; 
    15   - sudo apt-get install python-numpy python-scipy 
    16 #  - sudo apt-get install opencl-headers 
     37  - echo $TRAVIS_OS_NAME 
     38 
     39  - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then 
     40      wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; 
     41    fi; 
     42  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then 
     43      wget https://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O miniconda.sh; 
     44    fi; 
     45 
     46  - bash miniconda.sh -b -p $HOME/miniconda 
     47  - export PATH="$HOME/miniconda/bin:$PATH" 
     48  - hash -r 
     49  - conda update --yes conda 
     50 
     51  # Useful for debugging any issues with conda 
     52  - conda info -a 
     53 
     54  - conda install --yes python=$PY $NUMPYSPEC scipy cython mako cffi 
     55 
     56  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then 
     57      pip install pyopencl; 
     58    fi; 
    1759 
    1860install: 
    19 #  - pip install pyopencl 
    2061  - pip install bumps 
    2162  - pip install unittest-xml-reporting 
     63 
    2264script: 
    23   - export WORKSPACE=/home/travis/build/SasView/sasmodels/ 
    24   - python -m sasmodels.model_test dll all 
     65- export WORKSPACE=/home/travis/build/SasView/sasmodels/ 
     66- python -m sasmodels.model_test dll all 
     67 
     68notifications: 
     69  slack: 
     70    secure: xNAUeSu1/it/x9Q2CSg79aw1LLc7d6mLpcqSCTeKROp71RhkFf8VjJnJm/lEbKHNC8yj5H9UHrz5DmzwJzI+6oMt4NdEeS6WvGhwGY/wCt2IcJKxw0vj1DAU04qFMS041Khwclo6jIqm76DloinXvmvsS+K/nSyQkF7q4egSlwA= 

  • LICENSE.txt

    rc8de1bd rf68e2a5  
    1 This program is in the public domain. 
     1Copyright (c) 2009-2017, SasView Developers 
     2All rights reserved. 
    23 
    3 Individual files may be copyright different authors, with licences that 
    4 allow modification and redistribution in source or binary form with or 
    5 without modification, and a disclaimer of warranty. 
     4Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 
     5 
     61. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
     7 
     82. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 
     9 
     103. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. 
     11 
     12THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 

  • example/oriented_usans.py

    rea75043 r1cd24b4  
    1919    scale=0.08, background=0, 
    2020    sld=.291, sld_solvent=7.105, 
    21     r_polar=1800, r_polar_pd=0.222296, r_polar_pd_n=0, 
    22     r_equatorial=2600, r_equatorial_pd=0.28, r_equatorial_pd_n=0, 
     21    radius_polar=1800, radius_polar_pd=0.222296, radius_polar_pd_n=0, 
     22    radius_equatorial=2600, radius_equatorial_pd=0.28, radius_equatorial_pd_n=0, 
    2323    theta=60, theta_pd=0, theta_pd_n=0, 
    2424    phi=60, phi_pd=0, phi_pd_n=0, 
    … …  
    2626 
    2727# SET THE FITTING PARAMETERS 
    28 model.r_polar.range(1000, 10000) 
    29 model.r_equatorial.range(1000, 10000) 
     28model.radius_polar.range(1000, 10000) 
     29model.radius_equatorial.range(1000, 10000) 
    3030model.theta.range(0, 360) 
    3131model.phi.range(0, 360) 

  • explore/angular_pd.py

    r12eb36b r8267e0b  
    4747 
    4848def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    49     theta_center = radians(theta) 
     49    theta_center = radians(90-theta) 
    5050    phi_center = radians(phi) 
    5151    flow_center = radians(flow) 
    … …  
    137137                              radius=11., dist=dist) 
    138138        if not axis.startswith('d'): 
    139             ax.view_init(elev=theta, azim=phi) 
     139            ax.view_init(elev=90-theta if use_new else theta, azim=phi) 
    140140        plt.gcf().canvas.draw() 
    141141 

  • sasmodels/compare.py

    rf72d70a r630156b  
    3030 
    3131import sys 
     32import os 
    3233import math 
    3334import datetime 
    … …  
    4041from . import kerneldll 
    4142from . import exception 
    42 from .data import plot_theory, empty_data1D, empty_data2D 
     43from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4344from .direct_model import DirectModel 
    4445from .convert import revert_name, revert_pars, constrain_new_to_old 
    … …  
    7273    -1d*/-2d computes 1d or 2d data 
    7374    -preset*/-random[=seed] preset or random parameters 
    74     -mono/-poly* force monodisperse/polydisperse 
     75    -mono*/-poly force monodisperse or allow polydisperse demo parameters 
    7576    -magnetic/-nonmagnetic* suppress magnetism 
    7677    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    … …  
    8384    -edit starts the parameter explorer 
    8485    -default/-demo* use demo vs default parameters 
    85     -html shows the model docs instead of running the model 
     86    -help/-html shows the model docs instead of running the model 
    8687    -title="note" adds note to the plot title, after the model name 
     88    -data="path" uses q, dq from the data file 
    8789 
    8890Any two calculation engines can be selected for comparison: 
    … …  
    544546    'f32': '32', 
    545547    'f64': '64', 
     548    'float16': '16', 
     549    'float32': '32', 
     550    'float64': '64', 
     551    'float128': '128', 
    546552    'longdouble': '128', 
    547553} 
    … …  
    747753    comp = opts['engines'][1] if have_comp else None 
    748754    data = opts['data'] 
     755    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    749756 
    750757    # Plot if requested 
    751758    view = opts['view'] 
    752759    import matplotlib.pyplot as plt 
    753     if limits is None: 
     760    if limits is None and not use_data: 
    754761        vmin, vmax = np.Inf, -np.Inf 
    755762        if have_base: 
    … …  
    763770    if have_base: 
    764771        if have_comp: plt.subplot(131) 
    765         plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
     772        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    766773        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    767774        #cbar_title = "log I" 
    … …  
    769776        if have_base: plt.subplot(132) 
    770777        if not opts['is2d'] and have_base: 
    771             plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
    772         plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 
     778            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     779        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
    773780        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    774781        #cbar_title = "log I" 
    … …  
    784791            err[err>cutoff] = cutoff 
    785792        #err,errstr = base/comp,"ratio" 
    786         plot_theory(data, None, resid=err, view=errview, use_data=False) 
     793        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    787794        if view == 'linear': 
    788795            plt.xscale('linear') 
    … …  
    829836    'linear', 'log', 'q4', 
    830837    'hist', 'nohist', 
    831     'edit', 'html', 
     838    'edit', 'html', 'help', 
    832839    'demo', 'default', 
    833840    ]) 
    834841VALUE_OPTIONS = [ 
    835842    # Note: random is both a name option and a value option 
    836     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 
     843    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 
    837844    ] 
    838845 
    … …  
    940947        'cutoff'    : 0.0, 
    941948        'seed'      : -1,  # default to preset 
    942         'mono'      : False, 
     949        'mono'      : True, 
    943950        # Default to magnetic a magnetic moment is set on the command line 
    944951        'magnetic'  : False, 
    … …  
    951958        'html'      : False, 
    952959        'title'     : None, 
     960        'datafile'  : None, 
    953961    } 
    954962    engines = [] 
    … …  
    971979        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
    972980        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    973         elif arg.startswith('-title'):     opts['title'] = arg[7:] 
     981        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
     982        elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
    974983        elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
    975984        elif arg == '-preset':  opts['seed'] = -1 
    … …  
    9961005        elif arg == '-default':    opts['use_demo'] = False 
    9971006        elif arg == '-html':    opts['html'] = True 
     1007        elif arg == '-help':    opts['html'] = True 
    9981008    # pylint: enable=bad-whitespace 
    9991009 
    … …  
    11121122 
    11131123    # Create the computational engines 
    1114     data, _ = make_data(opts) 
     1124    if opts['datafile'] is not None: 
     1125        data = load_data(os.path.expanduser(opts['datafile'])) 
     1126    else: 
     1127        data, _ = make_data(opts) 
    11151128    if n1: 
    11161129        base = make_engine(model_info, data, engines[0], opts['cutoff']) 
    … …  
    11421155    show html docs for the model 
    11431156    """ 
    1144     import wx  # type: ignore 
    1145     from .generate import view_html_from_info 
    1146     app = wx.App() if wx.GetApp() is None else None 
    1147     view_html_from_info(opts['def'][0]) 
    1148     if app: app.MainLoop() 
    1149  
     1157    import os 
     1158    from .generate import make_html 
     1159    from . import rst2html 
     1160 
     1161    info = opts['def'][0] 
     1162    html = make_html(info) 
     1163    path = os.path.dirname(info.filename) 
     1164    url = "file://"+path.replace("\\","/")[2:]+"/" 
     1165    rst2html.view_html_qtapp(html, url) 
    11501166 
    11511167def explore(opts): 

  • sasmodels/core.py

    r5124c969 r650c6d2  
    250250        dtype = "longdouble" 
    251251    elif dtype == "half": 
    252         dtype = "f16" 
     252        dtype = "float16" 
    253253 
    254254    # Convert dtype string to numpy dtype. 

  • sasmodels/data.py

    r40a87fa r630156b  
    5151    from sas.sascalc.dataloader.loader import Loader  # type: ignore 
    5252    loader = Loader() 
    53     data = loader.load(filename) 
    54     if data is None: 
     53    # Allow for one part in multipart file 
     54    if '[' in filename: 
     55        filename, indexstr = filename[:-1].split('[') 
     56        index = int(indexstr) 
     57    else: 
     58        index = None 
     59    datasets = loader.load(filename) 
     60    if datasets is None: 
    5561        raise IOError("Data %r could not be loaded" % filename) 
     62    if not isinstance(datasets, list): 
     63        datasets = [datasets] 
     64    if index is None and len(datasets) > 1: 
     65        raise ValueError("Need to specify filename[index] for multipart data") 
     66    data = datasets[index if index is not None else 0] 
     67    if hasattr(data, 'x'): 
     68        data.qmin, data.qmax = data.x.min(), data.x.max() 
     69        data.mask = (np.isnan(data.y) if data.y is not None 
     70                     else np.zeros_like(data.x, dtype='bool')) 
     71    elif hasattr(data, 'qx_data'): 
     72        data.mask = ~data.mask 
    5673    return data 
    5774 
    … …  
    348365    # data, but they already handle the masking and graph markup already, so 
    349366    # do not repeat. 
    350     if hasattr(data, 'lam'): 
     367    if hasattr(data, 'isSesans') and data.isSesans: 
    351368        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    352369    elif hasattr(data, 'qx_data'): 
    … …  
    376393    *Iq_calc* is the raw theory values without resolution smearing 
    377394    """ 
    378     if hasattr(data, 'lam'): 
     395    if hasattr(data, 'isSesans') and data.isSesans: 
    379396        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    380397    elif hasattr(data, 'qx_data'): 
    … …  
    446463            if view is 'log': 
    447464                mtheory[mtheory <= 0] = masked 
    448             plt.plot(data.x, scale*mtheory, '-', hold=True) 
     465            plt.plot(data.x, scale*mtheory, '-') 
    449466            all_positive = all_positive and (mtheory > 0).all() 
    450467            some_present = some_present or (mtheory.count() > 0) 
    … …  
    453470            plt.ylim(*limits) 
    454471 
    455         plt.xscale('linear' if not some_present or non_positive_x  else view) 
     472        plt.xscale('linear' if not some_present or non_positive_x 
     473                   else view if view is not None 
     474                   else 'log') 
    456475        plt.yscale('linear' 
    457476                   if view == 'q4' or not some_present or not all_positive 
    458                    else view) 
     477                   else view if view is not None 
     478                   else 'log') 
    459479        plt.xlabel("$q$/A$^{-1}$") 
    460480        plt.ylabel('$I(q)$') 
     481        title = ("data and model" if use_theory and use_data 
     482                 else "data" if use_data 
     483                 else "model") 
     484        plt.title(title) 
    461485 
    462486    if use_calc: 
    … …  
    478502        if num_plots > 1: 
    479503            plt.subplot(1, num_plots, use_calc + 2) 
    480         plt.plot(data.x, mresid, '-') 
     504        plt.plot(data.x, mresid, '.') 
    481505        plt.xlabel("$q$/A$^{-1}$") 
    482506        plt.ylabel('residuals') 
    483         plt.xscale('linear' if not some_present or non_positive_x else view) 
     507        plt.xscale('linear') 
     508        plt.title('(model - Iq)/dIq') 
    484509 
    485510 
    … …  
    508533        if theory is not None: 
    509534            if is_tof: 
    510                 plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-', hold=True) 
     535                plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-') 
    511536            else: 
    512                 plt.plot(data.x, theory, '-', hold=True) 
     537                plt.plot(data.x, theory, '-') 
    513538        if limits is not None: 
    514539            plt.ylim(*limits) 

  • sasmodels/direct_model.py

    rb397165 ra769b54  
    192192 
    193193        # interpret data 
    194         if hasattr(data, 'lam'): 
     194        if hasattr(data, 'isSesans') and data.isSesans: 
    195195            self.data_type = 'sesans' 
    196196        elif hasattr(data, 'qx_data'): 

  • sasmodels/generate.py

    r1e7b0db0 rbb4b509  
    492492 
    493493def test_tag_float(): 
    494  
    495     cases=""" 
     494    """check that floating point constants are properly identified and tagged with 'f'""" 
     495 
     496    cases = """ 
    496497ZP  : 0. 
    497498ZPF : 0.0,0.01,0.1 
    … …  
    519520""" 
    520521 
    521     output=""" 
     522    output = """ 
    522523ZP  : 0.f 
    523524ZPF : 0.0f,0.01f,0.1f 
    … …  
    611612    # type: (str, List[Parameter]) -> List[str] 
    612613    """ 
    613     Return a list of *prefix.parameter* from parameter items. 
     614    Return a list of *prefix+parameter* from parameter items. 
     615 
     616    *prefix* should be "v." if v is a struct. 
    614617    """ 
    615618    return [p.as_call_reference(prefix) for p in pars] 
    … …  
    733736        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
    734737 
    735     magpars = [k-2 for k,p in enumerate(partable.call_parameters) 
     738    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    736739               if p.type == 'sld'] 
    737740 
    … …  
    742745    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    743746    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    744     for k,v in enumerate(magpars[:3]): 
     747    for k, v in enumerate(magpars[:3]): 
    745748        source.append("#define MAGNETIC_PAR%d %d"%(k+1, v)) 
    746749 
    … …  
    779782        "#undef CALL_IQ", 
    780783        "#undef KERNEL_NAME", 
    781          ] 
     784    ] 
    782785 
    783786    imagnetic = [ 
    … …  
    872875 
    873876# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog 
    874 RST_PROLOG = """\ 
     877RST_PROLOG = r"""\ 
    875878.. |Ang| unicode:: U+212B 
    876879.. |Ang^-1| replace:: |Ang|\ :sup:`-1` 
    … …  
    928931    from . import rst2html 
    929932    url = "file://"+dirname(info.filename)+"/" 
    930     rst2html.wxview(make_html(info), url=url) 
     933    rst2html.view_html(make_html(info), url=url) 
    931934 
    932935def demo_time(): 

  • sasmodels/kernel_header.c

    r1e7b0db0 rbb4b509  
    1414#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
    1515#  endif 
    16    // Intel CPU on Mac gives strange values for erf(); also on the tested 
     16   // Intel CPU on Mac gives strange values for erf(); on the verified 
    1717   // platforms (intel, nvidia, amd), the cephes erf() is significantly 
    1818   // faster than that available in the native OpenCL. 
    … …  
    5757         typedef int int32_t; 
    5858         #include <math.h> 
    59          // TODO: test isnan 
     59         // TODO: check isnan is correct 
    6060         inline double _isnan(double x) { return x != x; } // hope this doesn't optimize away! 
    6161         #undef isnan 
    … …  
    148148inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
     150// To rotate from the canonical position to theta, phi, psi, first rotate by 
     151// psi about the major axis, oriented along z, which is a rotation in the 
     152// detector plane xy. Next rotate by theta about the y axis, aligning the major 
     153// axis in the xz plane. Finally, rotate by phi in the detector plane xy. 
     154// To compute the scattering, undo these rotations in reverse order: 
     155//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 
     156// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 
     157// vector in the q direction. 
     158// To change between counterclockwise and clockwise rotation, change the 
     159// sign of phi and psi. 
     160 
    150161#if 1 
    151162//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
    … …  
    166177#endif 
    167178 
     179#if 1 
     180#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 
     181    q = sqrt(qx*qx + qy*qy); \ 
     182    const double qxhat = qx/q; \ 
     183    const double qyhat = qy/q; \ 
     184    double sin_theta, cos_theta; \ 
     185    double sin_phi, cos_phi; \ 
     186    double sin_psi, cos_psi; \ 
     187    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     188    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     189    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     190    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 
     191         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 
     192    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 
     193         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 
     194    zhat = qxhat*(-sin_theta*cos_phi) \ 
     195         + qyhat*(-sin_theta*sin_phi); \ 
     196    } while (0) 
     197#else 
     198// SasView 3.x definition of orientation 
    168199#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
    169200    q = sqrt(qx*qx + qy*qy); \ 
    … …  
    180211    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
    181212    } while (0) 
     213#endif 

  • sasmodels/kerneldll.py

    r886fa25 rb2f1d2f  
    9797 
    9898if "SAS_COMPILER" in os.environ: 
    99     compiler = os.environ["SAS_COMPILER"] 
     99    COMPILER = os.environ["SAS_COMPILER"] 
    100100elif os.name == 'nt': 
    101     # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment 
    102     # and we can use the MSVC compiler.  Otherwise, if tinycc is available 
    103     # the use it.  Otherwise, hope that mingw is available. 
    104     if "VCINSTALLDIR" in os.environ: 
    105         compiler = "msvc" 
    106     elif tinycc: 
    107         compiler = "tinycc" 
     101    if tinycc is not None: 
     102        COMPILER = "tinycc" 
     103    elif "VCINSTALLDIR" in os.environ: 
     104        # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment 
     105        # and we can use the MSVC compiler.  Otherwise, if tinycc is available 
     106        # the use it.  Otherwise, hope that mingw is available. 
     107        COMPILER = "msvc" 
    108108    else: 
    109         compiler = "mingw" 
     109        COMPILER = "mingw" 
    110110else: 
    111     compiler = "unix" 
     111    COMPILER = "unix" 
    112112 
    113113ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86 
    114 if compiler == "unix": 
     114if COMPILER == "unix": 
    115115    # Generic unix compile 
    116116    # On mac users will need the X code command line tools installed 
    … …  
    123123        """unix compiler command""" 
    124124        return CC + [source, "-o", output, "-lm"] 
    125 elif compiler == "msvc": 
     125elif COMPILER == "msvc": 
    126126    # Call vcvarsall.bat before compiling to set path, headers, libs, etc. 
    127127    # MSVC compiler is available, so use it.  OpenMP requires a copy of 
    … …  
    139139        """MSVC compiler command""" 
    140140        return CC + ["/Tp%s"%source] + LN + ["/OUT:%s"%output] 
    141 elif compiler == "tinycc": 
     141elif COMPILER == "tinycc": 
    142142    # TinyCC compiler. 
    143143    CC = [tinycc.TCC] + "-shared -rdynamic -Wall".split() 
    … …  
    145145        """tinycc compiler command""" 
    146146        return CC + [source, "-o", output] 
    147 elif compiler == "mingw": 
     147elif COMPILER == "mingw": 
    148148    # MinGW compiler. 
    149149    CC = "gcc -shared -std=c99 -O2 -Wall".split() 

  • sasmodels/kernelpy.py

    rb397165 rdbfd471  
    235235            # exclude all q for that NaN.  Even better would be to have an 
    236236            # INVALID expression like the C models, but that is too expensive. 
    237             Iq = form() 
     237            Iq = np.asarray(form(), 'd') 
    238238            if np.isnan(Iq).any(): continue 
    239239 

  • sasmodels/model_test.py

    r598857b rbb4b509  
    8585        skip = [] 
    8686    for model_name in models: 
    87         if model_name in skip: continue 
     87        if model_name in skip: 
     88            continue 
    8889        model_info = load_model_info(model_name) 
    8990 
    … …  
    239240            multiple = [test for test in self.info.tests 
    240241                        if isinstance(test[2], list) 
    241                             and not all(result is None for result in test[2])] 
     242                        and not all(result is None for result in test[2])] 
    242243            tests_has_1D_multiple = any(isinstance(test[1][0], float) 
    243244                                        for test in multiple) 
    … …  
    262263            user_pars, x, y = test 
    263264            pars = expand_pars(self.info.parameters, user_pars) 
     265            invalid = invalid_pars(self.info.parameters, pars) 
     266            if invalid: 
     267                raise ValueError("Unknown parameters in test: " + ", ".join(invalid)) 
    264268 
    265269            if not isinstance(y, list): 
    … …  
    305309 
    306310    return ModelTestCase 
     311 
     312def invalid_pars(partable, pars): 
     313    # type: (ParameterTable, Dict[str, float]) 
     314    """ 
     315    Return a list of parameter names that are not part of the model. 
     316    """ 
     317    names = set(p.id for p in partable.call_parameters) 
     318    invalid = [] 
     319    for par in sorted(pars.keys()): 
     320        parts = par.split('_pd') 
     321        if len(parts) > 1 and parts[1] not in ("", "_n", "nsigma", "type"): 
     322            invalid.append(par) 
     323            continue 
     324        if parts[0] not in names: 
     325            invalid.append(par) 
     326    return invalid 
     327 
    307328 
    308329def is_near(target, actual, digits=5): 

  • sasmodels/modelinfo.py

    rf88e248 rbb4b509  
    101101                limits = (float(low), float(high)) 
    102102            except Exception: 
    103                 raise ValueError("invalid limits for %s: %r"%(name,user_limits)) 
     103                raise ValueError("invalid limits for %s: %r"%(name, user_limits)) 
    104104            if low >= high: 
    105105                raise ValueError("require lower limit < upper limit") 
    … …  
    342342    def as_call_reference(self, prefix=""): 
    343343        # type: (str) -> str 
     344        """ 
     345        Return *prefix* + parameter name.  For struct references, use "v." 
     346        as the prefix. 
     347        """ 
    344348        # Note: if the parameter is a struct type, then we will need to use 
    345349        # &prefix+id.  For scalars and vectors we can just use prefix+id. 
    … …  
    420424        self.npars = sum(p.length for p in self.kernel_parameters) 
    421425        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
    422                              if p.type=='sld') 
     426                             if p.type == 'sld') 
    423427        self.nvalues = 2 + self.npars 
    424428        if self.nmagnetic: 
    … …  
    457461        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    458462                          for p in self.kernel_parameters) 
    459         self.magnetism_index = [k for k,p in enumerate(self.call_parameters) 
     463        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    460464                                if p.id.startswith('M0:')] 
    461465 
    … …  
    544548                              'magnetic', 'magnetic amplitude for '+p.description), 
    545549                    Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], 
    546                                'magnetic', 'magnetic latitude for '+p.description), 
     550                              'magnetic', 'magnetic latitude for '+p.description), 
    547551                    Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], 
    548                                'magnetic', 'magnetic longitude for '+p.description), 
     552                              'magnetic', 'magnetic longitude for '+p.description), 
    549553                ]) 
    550554        #print("call parameters", full_list) 
    … …  
    683687 
    684688    if (model_info.Iq is None 
    685         and model_info.Iqxy is None 
    686         and model_info.Imagnetic is None 
    687         and model_info.form_volume is None): 
     689            and model_info.Iqxy is None 
     690            and model_info.Imagnetic is None 
     691            and model_info.form_volume is None): 
    688692        return 
    689693 
    … …  
    843847    #: it to False because they require double precision calculations. 
    844848    single = None           # type: bool 
     849    #: True if the model can be run as an opencl model.  If for some reason 
     850    #: the model cannot be run in opencl (e.g., because the model passes 
     851    #: functions by reference), then set this to false. 
     852    opencl = None           # type: bool 
    845853    #: True if the model is a structure factor used to model the interaction 
    846854    #: between form factor models.  This will default to False if it is not 

  • sasmodels/models/barbell.py

    rfcb33e4 r9802ab3  
    6868The 2D scattering intensity is calculated similar to the 2D cylinder model. 
    6969 
    70 .. figure:: img/cylinder_angle_definition.jpg 
     70.. figure:: img/cylinder_angle_definition.png 
    7171 
    7272    Definition of the angles for oriented 2D barbells. 
    … …  
    8787* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    8888""" 
    89 from numpy import inf 
     89from numpy import inf, sin, cos, pi 
    9090 
    9191name = "barbell" 
    … …  
    108108              ["radius",      "Ang",         20, [0, inf],    "volume",      "Cylindrical bar radius"], 
    109109              ["length",      "Ang",        400, [0, inf],    "volume",      "Cylinder bar length"], 
    110               ["theta",       "degrees",     60, [-inf, inf], "orientation", "In plane angle"], 
    111               ["phi",         "degrees",     60, [-inf, inf], "orientation", "Out of plane angle"], 
     110              ["theta",       "degrees",     60, [-360, 360], "orientation", "Barbell axis to beam angle"], 
     111              ["phi",         "degrees",     60, [-360, 360], "orientation", "Rotation about beam"], 
    112112             ] 
    113113# pylint: enable=bad-whitespace, line-too-long 
    … …  
    125125            phi_pd=15, phi_pd_n=0, 
    126126           ) 
     127q = 0.1 
     128# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     129qx = q*cos(pi/6.0) 
     130qy = q*sin(pi/6.0) 
     131tests = [[{}, 0.075, 25.5691260532], 
     132        [{'theta':80., 'phi':10.}, (qx, qy), 3.04233067789], 
     133        ] 
     134del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/bcc_paracrystal.c

    r4962519 r50beefe  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, cos_a1, cos_a2, cos_a3; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
     92    double q, zhat, yhat, xhat; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    9494 
    95     const double a1 = +cos_a3 - cos_a1 + cos_a2; 
    96     const double a2 = +cos_a3 + cos_a1 - cos_a2; 
    97     const double a3 = -cos_a3 + cos_a1 + cos_a2; 
     95    const double a1 = +xhat - zhat + yhat; 
     96    const double a2 = +xhat + zhat - yhat; 
     97    const double a3 = -xhat + zhat + yhat; 
    9898 
    9999    const double qd = 0.5*q*dnn; 

  • sasmodels/models/bcc_paracrystal.py

    r925ad6e r69e1afc  
    7979be accurate. 
    8080 
    81 .. figure:: img/bcc_angle_definition.png 
     81.. figure:: img/parallelepiped_angle_definition.png 
    8282 
    83     Orientation of the crystal with respect to the scattering plane. 
     83    Orientation of the crystal with respect to the scattering plane, when  
     84    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8485 
    8586References 
    … …  
    99100""" 
    100101 
    101 from numpy import inf 
     102from numpy import inf, pi 
    102103 
    103104name = "bcc_paracrystal" 
    … …  
    122123              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "sld",         "Particle scattering length density"], 
    123124              ["sld_solvent", "1e-6/Ang^2",  1,    [-inf, inf], "sld",         "Solvent scattering length density"], 
    124               ["theta",       "degrees",    60,    [-inf, inf], "orientation", "In plane angle"], 
    125               ["phi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"], 
    126               ["psi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"] 
     125              ["theta",       "degrees",    60,    [-360, 360], "orientation", "c axis to beam angle"], 
     126              ["phi",         "degrees",    60,    [-360, 360], "orientation", "rotation about beam"], 
     127              ["psi",         "degrees",    60,    [-360, 360], "orientation", "rotation about c axis"] 
    127128             ] 
    128129# pylint: enable=bad-whitespace, line-too-long 
    … …  
    141142    psi_pd=15, psi_pd_n=0, 
    142143    ) 
     144# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     145# add 2d test later 
     146q =4.*pi/220. 
     147tests = [ 
     148    [{ }, 
     149     [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]], 
     150    [{'theta':20.0,'phi':30,'psi':40.0},(-0.017,0.035),2082.20264399 ], 
     151    [{'theta':20.0,'phi':30,'psi':40.0},(-0.081,0.011),0.436323144781 ] 
     152    ] 

  • sasmodels/models/capped_cylinder.py

    rfcb33e4 r9802ab3  
    7171The 2D scattering intensity is calculated similar to the 2D cylinder model. 
    7272 
    73 .. figure:: img/cylinder_angle_definition.jpg 
     73.. figure:: img/cylinder_angle_definition.png 
    7474 
    7575    Definition of the angles for oriented 2D cylinders. 
    … …  
    9191 
    9292""" 
    93 from numpy import inf 
     93from numpy import inf, sin, cos, pi 
    9494 
    9595name = "capped_cylinder" 
    … …  
    129129              ["radius_cap", "Ang",     20, [0, inf],    "volume", "Cap radius"], 
    130130              ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"], 
    131               ["theta",      "degrees", 60, [-inf, inf], "orientation", "inclination angle"], 
    132               ["phi",        "degrees", 60, [-inf, inf], "orientation", "deflection angle"], 
     131              ["theta",      "degrees", 60, [-360, 360], "orientation", "cylinder axis to beam angle"], 
     132              ["phi",        "degrees", 60, [-360, 360], "orientation", "rotation about beam"], 
    133133             ] 
    134134# pylint: enable=bad-whitespace, line-too-long 
    … …  
    145145            theta_pd=15, theta_pd_n=45, 
    146146            phi_pd=15, phi_pd_n=1) 
     147q = 0.1 
     148# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     149qx = q*cos(pi/6.0) 
     150qy = q*sin(pi/6.0) 
     151tests = [[{}, 0.075, 26.0698570695], 
     152        [{'theta':80., 'phi':10.}, (qx, qy), 0.561811990502], 
     153        ] 
     154del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/core_shell_bicelle.c

    r592343f rb260926  
    3030 
    3131static double 
    32 bicelle_kernel(double qq, 
     32bicelle_kernel(double q, 
    3333              double rad, 
    3434              double radthick, 
    3535              double facthick, 
    36               double length, 
     36              double halflength, 
    3737              double rhoc, 
    3838              double rhoh, 
    … …  
    4242              double cos_alpha) 
    4343{ 
    44     double si1,si2,be1,be2; 
    45  
    4644    const double dr1 = rhoc-rhoh; 
    4745    const double dr2 = rhor-rhosolv; 
    4846    const double dr3 = rhoh-rhor; 
    49     const double vol1 = M_PI*rad*rad*(2.0*length); 
    50     const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 
    51     const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 
    52     double besarg1 = qq*rad*sin_alpha; 
    53     double besarg2 = qq*(rad+radthick)*sin_alpha; 
    54     double sinarg1 = qq*length*cos_alpha; 
    55     double sinarg2 = qq*(length+facthick)*cos_alpha; 
     47    const double vol1 = M_PI*square(rad)*2.0*(halflength); 
     48    const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 
     49    const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 
    5650 
    57     be1 = sas_2J1x_x(besarg1); 
    58     be2 = sas_2J1x_x(besarg2); 
    59     si1 = sas_sinx_x(sinarg1); 
    60     si2 = sas_sinx_x(sinarg2); 
     51    const double be1 = sas_2J1x_x(q*(rad)*sin_alpha); 
     52    const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha); 
     53    const double si1 = sas_sinx_x(q*(halflength)*cos_alpha); 
     54    const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha); 
    6155 
    6256    const double t = vol1*dr1*si1*be1 + 
    … …  
    6458                     vol3*dr3*si2*be1; 
    6559 
    66     const double retval = t*t*sin_alpha; 
     60    const double retval = t*t; 
    6761 
    6862    return retval; 
    … …  
    7165 
    7266static double 
    73 bicelle_integration(double qq, 
     67bicelle_integration(double q, 
    7468                   double rad, 
    7569                   double radthick, 
    … …  
    8377    // set up the integration end points 
    8478    const double uplim = M_PI_4; 
    85     const double halfheight = 0.5*length; 
     79    const double halflength = 0.5*length; 
    8680 
    8781    double summ = 0.0; 
    … …  
    9084        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    9185        SINCOS(alpha, sin_alpha, cos_alpha); 
    92         double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
    93                              halfheight, rhoc, rhoh, rhor, rhosolv, 
     86        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
     87                             halflength, rhoc, rhoh, rhor, rhosolv, 
    9488                             sin_alpha, cos_alpha); 
    95         summ += yyy; 
     89        summ += yyy*sin_alpha; 
    9690    } 
    9791 
    … …  
    119113    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
    120114                           0.5*length, core_sld, face_sld, rim_sld, 
    121                            solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 
    122  
    123     answer *= 1.0e-4; 
    124  
    125     return answer; 
     115                           solvent_sld, sin_alpha, cos_alpha); 
     116    return 1.0e-4*answer; 
    126117} 
    127118 

  • sasmodels/models/core_shell_bicelle.py

    r8afefae r9802ab3  
    4747.. math:: 
    4848 
    49         \begin{align}     
     49    \begin{align}     
    5050    F(Q,\alpha) = &\bigg[  
    5151    (\rho_c - \rho_f) V_c \frac{2J_1(QRsin \alpha)}{QRsin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\ 
    … …  
    6363cylinders is then given by integrating over all possible $\theta$ and $\phi$. 
    6464 
    65 The *theta* and *phi* parameters are not used for the 1D output. 
     65For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data,  
     66see the :ref:`cylinder` model for further information. 
    6667Our implementation of the scattering kernel and the 1D scattering intensity 
    6768use the c-library from NIST. 
    6869 
    69 .. figure:: img/cylinder_angle_definition.jpg 
     70.. figure:: img/cylinder_angle_definition.png 
    7071 
    71     Definition of the angles for the oriented core shell bicelle tmodel. 
     72    Definition of the angles for the oriented core shell bicelle model, 
     73    note that the cylinder axis of the bicelle starts along the beam direction 
     74    when $\theta  = \phi = 0$. 
    7275 
    7376 
    … …  
    8891""" 
    8992 
    90 from numpy import inf, sin, cos 
     93from numpy import inf, sin, cos, pi 
    9194 
    9295name = "core_shell_bicelle" 
    … …  
    135138    ["sld_rim",        "1e-6/Ang^2", 4, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    136139    ["sld_solvent",    "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Solvent scattering length density"], 
    137     ["theta",          "degrees",   90, [-inf, inf], "orientation", "In plane angle"], 
    138     ["phi",            "degrees",    0, [-inf, inf], "orientation", "Out of plane angle"], 
     140    ["theta",          "degrees",   90, [-360, 360], "orientation", "cylinder axis to beam angle"], 
     141    ["phi",            "degrees",    0, [-360, 360], "orientation", "rotation about beam"] 
    139142    ] 
    140143 
    … …  
    155158            theta=90, 
    156159            phi=0) 
     160q = 0.1 
     161# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     162qx = q*cos(pi/6.0) 
     163qy = q*sin(pi/6.0) 
     164tests = [[{}, 0.05, 7.4883545957], 
     165        [{'theta':80., 'phi':10.}, (qx, qy), 2.81048892474 ] 
     166        ] 
     167del qx, qy  # not necessary to delete, but cleaner 
    157168 
    158 #qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r592343f rdedcf34  
    1 double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length); 
    2 double Iq(double q, 
    3           double radius, 
    4           double x_core, 
    5           double thick_rim, 
    6           double thick_face, 
    7           double length, 
    8           double core_sld, 
    9           double face_sld, 
    10           double rim_sld, 
    11           double solvent_sld); 
    12  
    13  
    14 double Iqxy(double qx, double qy, 
    15           double radius, 
    16           double x_core, 
    17           double thick_rim, 
    18           double thick_face, 
    19           double length, 
    20           double core_sld, 
    21           double face_sld, 
    22           double rim_sld, 
    23           double solvent_sld, 
    24           double theta, 
    25           double phi, 
    26           double psi); 
    27  
    281// NOTE that "length" here is the full height of the core! 
    29 double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length) 
     2static double 
     3form_volume(double r_minor, 
     4        double x_core, 
     5        double thick_rim, 
     6        double thick_face, 
     7        double length) 
    308{ 
    31     return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face); 
     9    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 
    3210} 
    3311 
    34 double  
    35                 Iq(double qq, 
    36                    double rad, 
    37                    double x_core, 
    38                    double radthick, 
    39                    double facthick, 
    40                    double length, 
    41                    double rhoc, 
    42                    double rhoh, 
    43                    double rhor, 
    44                    double rhosolv) 
     12static double 
     13Iq(double q, 
     14        double r_minor, 
     15        double x_core, 
     16        double thick_rim, 
     17        double thick_face, 
     18        double length, 
     19        double rhoc, 
     20        double rhoh, 
     21        double rhor, 
     22        double rhosolv) 
    4523{ 
    4624    double si1,si2,be1,be2; 
    4725     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    48      // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle 
     26     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
    4927     //    const double uplim = M_PI_4; 
    5028    const double halfheight = 0.5*length; 
    … …  
    5533    //const double vbj=M_PI; 
    5634 
    57     const double radius_major = rad * x_core; 
    58     const double rA = 0.5*(square(radius_major) + square(rad)); 
    59     const double rB = 0.5*(square(radius_major) - square(rad)); 
    60     const double dr1 = (rhoc-rhoh)   *M_PI*rad*radius_major*(2.0*halfheight);; 
    61     const double dr2 = (rhor-rhosolv)*M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick); 
    62     const double dr3 = (rhoh-rhor)   *M_PI*rad*radius_major*2.0*(halfheight+facthick); 
    63     //const double vol1 = M_PI*rad*radius_major*(2.0*halfheight); 
    64     //const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick); 
    65     //const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 
     35    const double r_major = r_minor * x_core; 
     36    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
     37    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
     38    const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);; 
     39    const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     40    const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     41    //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     42    //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     43    //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    6644 
    6745    //initialize integral 
    … …  
    7452        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    7553        double inner_sum=0; 
    76         double sinarg1 = qq*halfheight*cos_alpha; 
    77         double sinarg2 = qq*(halfheight+facthick)*cos_alpha; 
     54        double sinarg1 = q*halfheight*cos_alpha; 
     55        double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    7856        si1 = sas_sinx_x(sinarg1); 
    7957        si2 = sas_sinx_x(sinarg2); 
    … …  
    8260            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    8361            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    84             const double rr = sqrt(rA - rB*cos(beta)); 
    85             double besarg1 = qq*rr*sin_alpha; 
    86             double besarg2 = qq*(rr+radthick)*sin_alpha; 
     62            const double rr = sqrt(r2A - r2B*cos(beta)); 
     63            double besarg1 = q*rr*sin_alpha; 
     64            double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    8765            be1 = sas_2J1x_x(besarg1); 
    8866            be2 = sas_2J1x_x(besarg2); 
    … …  
    9876} 
    9977 
    100 double  
     78static double 
    10179Iqxy(double qx, double qy, 
    102           double rad, 
     80          double r_minor, 
    10381          double x_core, 
    104           double radthick, 
    105           double facthick, 
     82          double thick_rim, 
     83          double thick_face, 
    10684          double length, 
    10785          double rhoc, 
    … …  
    11492{ 
    11593       // THIS NEEDS TESTING 
    116     double qq, cos_val, cos_mu, cos_nu; 
    117     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu); 
     94    double q, xhat, yhat, zhat; 
     95    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    11896    const double dr1 = rhoc-rhoh; 
    11997    const double dr2 = rhor-rhosolv; 
    12098    const double dr3 = rhoh-rhor; 
    121     const double radius_major = rad*x_core; 
     99    const double r_major = r_minor*x_core; 
    122100    const double halfheight = 0.5*length; 
    123     const double vol1 = M_PI*rad*radius_major*length; 
    124     const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick); 
    125     const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 
     101    const double vol1 = M_PI*r_minor*r_major*length; 
     102    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     103    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    126104 
    127     // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    128     // Given:    radius_major = r_ratio * radius_minor   
    129     // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle 
    130     const double r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu); 
    131     const double be1 = sas_2J1x_x( qq*r ); 
    132     const double be2 = sas_2J1x_x( qq*(r + radthick ) ); 
    133     const double si1 = sas_sinx_x( qq*halfheight*cos_val ); 
    134     const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val ); 
     105    // Compute effective radius in rotated coordinates 
     106    const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat)); 
     107    const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat) 
     108                                   + square((r_minor+thick_rim)*yhat)); 
     109    const double be1 = sas_2J1x_x( q*r_hat ); 
     110    const double be2 = sas_2J1x_x( q*rshell_hat ); 
     111    const double si1 = sas_sinx_x( q*halfheight*zhat ); 
     112    const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat ); 
    135113    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    136     //const double vol = form_volume(radius_minor, r_ratio, length); 
    137114    return 1.0e-4 * Aq; 
    138115} 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r8afefae r9802ab3  
    7676bicelles is then given by integrating over all possible $\alpha$ and $\psi$. 
    7777 
    78 For oriented bicellles the *theta*, *phi* and *psi* orientation parameters only appear when fitting 2D data,  
     78For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
    7979see the :ref:`elliptical-cylinder` model for further information. 
    8080 
    8181 
    82 .. figure:: img/elliptical_cylinder_angle_definition.jpg 
     82.. figure:: img/elliptical_cylinder_angle_definition.png 
    8383 
    84     Definition of the angles for the oriented core_shell_bicelle_elliptical model. 
    85     Note that *theta* and *phi* are currently defined differently to those for the core_shell_bicelle model. 
     84    Definition of the angles for the oriented core_shell_bicelle_elliptical particles.    
     85 
    8686 
    8787 
    … …  
    9999""" 
    100100 
    101 from numpy import inf, sin, cos 
     101from numpy import inf, sin, cos, pi 
    102102 
    103103name = "core_shell_bicelle_elliptical" 
    … …  
    119119    ["radius",         "Ang",       30, [0, inf],    "volume",      "Cylinder core radius"], 
    120120    ["x_core",        "None",       3,  [0, inf],    "volume",      "axial ratio of core, X = r_polar/r_equatorial"], 
    121     ["thick_rim",  "Ang",        8, [0, inf],    "volume",      "Rim shell thickness"], 
    122     ["thick_face", "Ang",       14, [0, inf],    "volume",      "Cylinder face thickness"], 
    123     ["length",         "Ang",      50, [0, inf],    "volume",      "Cylinder length"], 
     121    ["thick_rim",  "Ang",            8, [0, inf],    "volume",      "Rim shell thickness"], 
     122    ["thick_face", "Ang",           14, [0, inf],    "volume",      "Cylinder face thickness"], 
     123    ["length",         "Ang",       50, [0, inf],    "volume",      "Cylinder length"], 
    124124    ["sld_core",       "1e-6/Ang^2", 4, [-inf, inf], "sld",         "Cylinder core scattering length density"], 
    125125    ["sld_face",       "1e-6/Ang^2", 7, [-inf, inf], "sld",         "Cylinder face scattering length density"], 
    126126    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    127127    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
    128     ["theta",          "degrees",   90, [-360, 360], "orientation", "In plane angle"], 
    129     ["phi",            "degrees",    0, [-360, 360], "orientation", "Out of plane angle"], 
    130     ["psi",            "degrees",    0, [-360, 360], "orientation", "Major axis angle relative to Q"], 
     128    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
     129    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
     130    ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"] 
    131131    ] 
    132132 
    … …  
    150150            psi=0) 
    151151 
    152 #qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 
     152q = 0.1 
     153# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
     154qx = q*cos(pi/6.0) 
     155qy = q*sin(pi/6.0) 
    153156 
    154157tests = [ 
    … …  
    159162    'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 
    160163    0.015, 286.540286], 
    161 ] 
     164#    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
     165        ] 
     166 
     167del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/core_shell_cylinder.py

    r8e68ea0 r9b79f29  
    7373""" 
    7474 
    75 from numpy import pi, inf 
     75from numpy import pi, inf, sin, cos 
    7676 
    7777name = "core_shell_cylinder" 
    … …  
    117117              ["length", "Ang", 400, [0, inf], "volume", 
    118118               "Cylinder length"], 
    119               ["theta", "degrees", 60, [-inf, inf], "orientation", 
    120                "In plane angle"], 
    121               ["phi", "degrees", 60, [-inf, inf], "orientation", 
    122                "Out of plane angle"], 
     119              ["theta", "degrees", 60, [-360, 360], "orientation", 
     120               "cylinder axis to beam angle"], 
     121              ["phi", "degrees",   60, [-360, 360], "orientation", 
     122               "rotation about beam"], 
    123123             ] 
    124124 
    … …  
    151151            theta_pd=15, theta_pd_n=45, 
    152152            phi_pd=15, phi_pd_n=1) 
    153  
     153q = 0.1 
     154# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     155qx = q*cos(pi/6.0) 
     156qy = q*sin(pi/6.0) 
     157tests = [[{}, 0.075, 10.8552692237], 
     158        [{}, (qx, qy), 0.444618752741 ], 
     159        ] 
     160del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/core_shell_ellipsoid.py

    r8e68ea0 r9802ab3  
    7777   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    7878 
     79For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     80see the :ref:`elliptical-cylinder` model for further information. 
    7981 
    8082References 
    … …  
    132134    ["sld_shell",     "1e-6/Ang^2", 1,   [-inf, inf], "sld",         "Shell scattering length density"], 
    133135    ["sld_solvent",   "1e-6/Ang^2", 6.3, [-inf, inf], "sld",         "Solvent scattering length density"], 
    134     ["theta",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"], 
    135     ["phi",           "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"], 
     136    ["theta",         "degrees",    0,   [-360, 360], "orientation", "elipsoid axis to beam angle"], 
     137    ["phi",           "degrees",    0,   [-360, 360], "orientation", "rotation about beam"], 
    136138    ] 
    137139# pylint: enable=bad-whitespace, line-too-long 

  • sasmodels/models/core_shell_parallelepiped.c

    r1e7b0db0 r92dfe0c  
    134134    double psi) 
    135135{ 
    136     double q, cos_val_a, cos_val_b, cos_val_c; 
    137     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
     136    double q, zhat, yhat, xhat; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    138138 
    139139    // cspkernel in csparallelepiped recoded here 
    … …  
    160160    double tc = length_a + 2.0*thick_rim_c; 
    161161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*cos_val_a); 
    163     double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 
    164     double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 
    165     double siAt = sas_sinx_x(0.5*q*ta*cos_val_a); 
    166     double siBt = sas_sinx_x(0.5*q*tb*cos_val_b); 
    167     double siCt = sas_sinx_x(0.5*q*tc*cos_val_c); 
     162    double siA = sas_sinx_x(0.5*q*length_a*xhat); 
     163    double siB = sas_sinx_x(0.5*q*length_b*yhat); 
     164    double siC = sas_sinx_x(0.5*q*length_c*zhat); 
     165    double siAt = sas_sinx_x(0.5*q*ta*xhat); 
     166    double siBt = sas_sinx_x(0.5*q*tb*yhat); 
     167    double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168168     
    169169 

  • sasmodels/models/core_shell_parallelepiped.py

    rcb0dc22 r9b79f29  
    66The thickness and the scattering length density of the shell or  
    77"rim" can be different on each (pair) of faces. However at this time 
    8 the model does **NOT** actually calculate a c face rim despite the presence of 
    9 the parameter. 
     8the 1D calculation does **NOT** actually calculate a c face rim despite the presence of 
     9the parameter. Some other aspects of the 1D calculation may be wrong. 
    1010 
    1111.. note:: 
    12    This model was originally ported from NIST IGOR macros. However,t is not 
    13    yet fully understood by the SasView developers and is currently review. 
     12   This model was originally ported from NIST IGOR macros. However, it is not 
     13   yet fully understood by the SasView developers and is currently under review. 
    1414 
    1515The form factor is normalized by the particle volume $V$ such that 
    … …  
    4848amplitudes of the core and shell, in the same manner as a core-shell model. 
    4949 
     50.. math:: 
     51 
     52    F_{a}(Q,\alpha,\beta)= 
     53    \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)}{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
     54    - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
     55    \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
     56    \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
    5057 
    5158.. note:: 
    5259 
     60    Why does t_B not appear in the above equation? 
    5361    For the calculation of the form factor to be valid, the sides of the solid 
    54     MUST be chosen such that** $A < B < C$. 
     62    MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
    5563    If this inequality is not satisfied, the model will not report an error, 
    5664    but the calculation will not be correct and thus the result wrong. 
    5765 
    5866FITTING NOTES 
    59 If the scale is set equal to the particle volume fraction, |phi|, the returned 
     67If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    6068value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    6169However, **no interparticle interference effects are included in this 
    … …  
    7381NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    7482based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    75 and length $(C+2t_C)$ values, and used as the effective radius 
    76 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     83and length $(C+2t_C)$ values, after appropriately 
     84sorting the three dimensions to give an oblate or prolate particle, to give an  
     85effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    7786 
    7887To provide easy access to the orientation of the parallelepiped, we define the 
    … …  
    8392*x*-axis of the detector. 
    8493 
    85 .. figure:: img/parallelepiped_angle_definition.jpg 
     94.. figure:: img/parallelepiped_angle_definition.png 
    8695 
    8796    Definition of the angles for oriented core-shell parallelepipeds. 
    8897 
    89 .. figure:: img/parallelepiped_angle_projection.jpg 
     98.. figure:: img/parallelepiped_angle_projection.png 
    9099 
    91100    Examples of the angles for oriented core-shell parallelepipeds against the 
    … …  
    112121 
    113122import numpy as np 
    114 from numpy import pi, inf, sqrt 
     123from numpy import pi, inf, sqrt, cos, sin 
    115124 
    116125name = "core_shell_parallelepiped" 
    … …  
    144153              ["thick_rim_c", "Ang", 10, [0, inf], "volume", 
    145154               "Thickness of C rim"], 
    146               ["theta", "degrees", 0, [-inf, inf], "orientation", 
    147                "In plane angle"], 
    148               ["phi", "degrees", 0, [-inf, inf], "orientation", 
    149                "Out of plane angle"], 
    150               ["psi", "degrees", 0, [-inf, inf], "orientation", 
    151                "Rotation angle around its own c axis against q plane"], 
     155              ["theta", "degrees", 0, [-360, 360], "orientation", 
     156               "c axis to beam angle"], 
     157              ["phi", "degrees", 0, [-360, 360], "orientation", 
     158               "rotation about beam"], 
     159              ["psi", "degrees", 0, [-360, 360], "orientation", 
     160               "rotation about c axis"], 
    152161             ] 
    153162 
    … …  
    186195            psi_pd=10, psi_pd_n=1) 
    187196 
    188 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     197# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     198qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    189199tests = [[{}, 0.2, 0.533149288477], 
    190200         [{}, [0.2], [0.533149288477]], 
    191          [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.032102135569], 
    192          [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.032102135569]], 
     201         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
     202         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
    193203        ] 
    194204del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/cylinder.py

    rb7e8b94 r9802ab3  
    3838 
    3939 
    40 Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with  
    41 $sin(\alpha)=\sqrt{1-u^2}$.  
     40Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with 
     41$sin(\alpha)=\sqrt{1-u^2}$. 
    4242 
    4343The output of the 1D scattering intensity function for randomly oriented 
    … …  
    6161.. _cylinder-angle-definition: 
    6262 
    63 .. figure:: img/cylinder_angle_definition.jpg 
     63.. figure:: img/cylinder_angle_definition.png 
    6464 
    65     Definition of the angles for oriented cylinders. 
     65    Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative  
     66    to the beam line coordinates, plus an indication of their orientation distributions  
     67    which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$  
     68    in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes. 
    6669 
    67 The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
     70.. figure:: img/cylinder_angle_projection.png 
     71 
     72    Examples for oriented cylinders. 
     73 
     74The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.  
     75On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
     76appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel  
     77to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully. 
     78(Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such  
     79situations, with very broad distributions, should still be approached with care.)  
    6880 
    6981Validation 
    … …  
    123135              ["length", "Ang", 400, [0, inf], "volume", 
    124136               "Cylinder length"], 
    125               ["theta", "degrees", 60, [-inf, inf], "orientation", 
    126                "latitude"], 
    127               ["phi", "degrees", 60, [-inf, inf], "orientation", 
    128                "longitude"], 
     137              ["theta", "degrees", 60, [-360, 360], "orientation", 
     138               "cylinder axis to beam angle"], 
     139              ["phi", "degrees",   60, [-360, 360], "orientation", 
     140               "rotation about beam"], 
    129141             ] 
    130142 
    … …  
    152164tests = [[{}, 0.2, 0.042761386790780453], 
    153165        [{}, [0.2], [0.042761386790780453]], 
    154 #  new coords     
     166#  new coords 
    155167        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
    156168        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 

  • sasmodels/models/ellipsoid.c

    r130d4c7 r3b571ae  
    33double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 static double 
    7 _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
    8 { 
    9     double ratio = radius_polar/radius_equatorial; 
    10     // Using ratio v = Rp/Re, we can expand the following to match the 
    11     // form given in Guinier (1955) 
    12     //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1)) 
    13     //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) ) 
    14     //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) ) 
    15     //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) ) 
    16     // 
    17     // Instead of using pythagoras we could pass in sin and cos; this may be 
    18     // slightly better for 2D which has already computed it, but it introduces 
    19     // an extra sqrt and square for 1-D not required by the current form, so 
    20     // leave it as is. 
    21     const double r = radius_equatorial 
    22                      * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
    23     const double f = sas_3j1x_x(q*r); 
    24  
    25     return f*f; 
    26 } 
    275 
    286double form_volume(double radius_polar, double radius_equatorial) 
    … …  
    3715    double radius_equatorial) 
    3816{ 
     17    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
     18    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
     19    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
     20    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
     21    // u-substitution of 
     22    //     u = sin, du = cos dT 
     23    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
     24    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
     25 
    3926    // translate a point in [-1,1] to a point in [0, 1] 
     27    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    4028    const double zm = 0.5; 
    4129    const double zb = 0.5; 
    4230    double total = 0.0; 
    4331    for (int i=0;i<76;i++) { 
    44         //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    45         const double cos_alpha = Gauss76Z[i]*zm + zb; 
    46         total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
     32        const double u = Gauss76Z[i]*zm + zb; 
     33        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
     34        const double f = sas_3j1x_x(q*r); 
     35        total += Gauss76Wt[i] * f * f; 
    4736    } 
    4837    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6251    double q, sin_alpha, cos_alpha; 
    6352    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    64     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
     53    const double r = sqrt(square(radius_equatorial*sin_alpha) 
     54                          + square(radius_polar*cos_alpha)); 
     55    const double f = sas_3j1x_x(q*r); 
    6556    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6657 
    67     return 1.0e-4 * form * s * s; 
     58    return 1.0e-4 * square(f * s); 
    6859} 
    6960 

  • sasmodels/models/ellipsoid.py

    r925ad6e r9b79f29  
    1818.. math:: 
    1919 
    20     F(q,\alpha) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)] 
    21                 - \cos[qr(R_p,R_e,\alpha)])} 
    22                 {[qr(R_p,R_e,\alpha)]^3} 
     20    F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3} 
    2321 
    24 and 
     22for 
    2523 
    2624.. math:: 
    2725 
    28     r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha 
    29         + R_p^2 \cos^2 \alpha \right]^{1/2} 
     26    r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2} 
    3027 
    3128 
    3229$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$, 
    33 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the 
    34 rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular 
    35 to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the 
    36 scattering length density difference between the scatterer and the solvent. 
     30$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar 
     31radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 
     32radius perpendicular to the rotational axis of the ellipsoid and 
     33$\Delta \rho$ (contrast) is the scattering length density difference between 
     34the scatterer and the solvent. 
    3735 
    38 For randomly oriented particles: 
     36For randomly oriented particles use the orientational average, 
    3937 
    4038.. math:: 
    4139 
    42    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
     40   \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha} 
    4341 
     42 
     43computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as 
     44 
     45.. math:: 
     46 
     47    \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du} 
     48 
     49with 
     50 
     51.. math:: 
     52 
     53    r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2} 
    4454 
    4555To provide easy access to the orientation of the ellipsoid, we define 
    … …  
    4858:ref:`cylinder orientation figure <cylinder-angle-definition>`. 
    4959For the ellipsoid, $\theta$ is the angle between the rotational axis 
    50 and the $z$ -axis. 
     60and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 
     61in the $xy$ plane. 
    5162 
    5263NB: The 2nd virial coefficient of the solid ellipsoid is calculated based 
    … …  
    90101than 500. 
    91102 
     103Model was also tested against the triaxial ellipsoid model with equal major 
     104and minor equatorial radii.  It is also consistent with the cyclinder model 
     105with polar radius equal to length and equatorial radius equal to radius. 
     106 
    92107References 
    93108---------- 
    … …  
    96111*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
    97112Plenum Press, New York, 1987. 
     113 
     114Authorship and Verification 
     115---------------------------- 
     116 
     117* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     118* **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014 
     119* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 
    98120""" 
    99121 
    100 from numpy import inf 
     122from numpy import inf, sin, cos, pi 
    101123 
    102124name = "ellipsoid" 
    … …  
    129151              ["radius_equatorial", "Ang", 400, [0, inf], "volume", 
    130152               "Equatorial radius"], 
    131               ["theta", "degrees", 60, [-inf, inf], "orientation", 
    132                "In plane angle"], 
    133               ["phi", "degrees", 60, [-inf, inf], "orientation", 
    134                "Out of plane angle"], 
     153              ["theta", "degrees", 60, [-360, 360], "orientation", 
     154               "ellipsoid axis to beam angle"], 
     155              ["phi", "degrees", 60, [-360, 360], "orientation", 
     156               "rotation about beam"], 
    135157             ] 
    136158 
    … …  
    139161def ER(radius_polar, radius_equatorial): 
    140162    import numpy as np 
    141  
     163# see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    142164    ee = np.empty_like(radius_polar) 
    143165    idx = radius_polar > radius_equatorial 
    … …  
    168190            theta_pd=15, theta_pd_n=45, 
    169191            phi_pd=15, phi_pd_n=1) 
     192q = 0.1 
     193# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     194qx = q*cos(pi/6.0) 
     195qy = q*sin(pi/6.0) 
     196tests = [[{}, 0.05, 54.8525847025], 
     197        [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 ], 
     198        ] 
     199del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/elliptical_cylinder.c

    r592343f r61104c8  
    6767     double theta, double phi, double psi) 
    6868{ 
    69     double q, cos_val, cos_mu, cos_nu; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 
     69    double q, xhat, yhat, zhat; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    7171 
    7272    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7373    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 
     74    const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    7575    const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*0.5*length*cos_val); 
     76    const double si = sas_sinx_x(q*zhat*0.5*length); 
    7777    const double Aq = be * si; 
    7878    const double delrho = sld - solvent_sld; 

  • sasmodels/models/elliptical_cylinder.py

    rfcb33e4 r9802ab3  
    5757define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$ 
    5858(see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle 
    59 $\Psi$ is the rotational angle around its own long_c axis against the $q$ plane. 
    60 For example, $\Psi = 0$ when the $r_\text{minor}$ axis is parallel to the 
    61 $x$ axis of the detector. 
     59$\Psi$ is the rotational angle around its own long_c axis.  
    6260 
    6361All angle parameters are valid and given only for 2D calculation; ie, an 
    6462oriented system. 
    6563 
    66 .. figure:: img/elliptical_cylinder_angle_definition.jpg 
     64.. figure:: img/elliptical_cylinder_angle_definition.png 
    6765 
    68     Definition of angles for 2D 
     66    Definition of angles for oriented elliptical cylinder, where axis_ratio is drawn >1, 
     67    and angle $\Psi$ is now a rotation around the axis of the cylinder. 
    6968 
    70 .. figure:: img/cylinder_angle_projection.jpg 
     69.. figure:: img/elliptical_cylinder_angle_projection.png 
    7170 
    7271    Examples of the angles for oriented elliptical cylinders against the 
    73     detector plane. 
     72    detector plane, with $\Psi$ = 0. 
     73 
     74The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.  
     75On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
     76appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the  
     77cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)  
     78The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to  
     79understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation  
     80distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)  
    7481 
    7582NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    … …  
    108115""" 
    109116 
    110 from numpy import pi, inf, sqrt 
     117from numpy import pi, inf, sqrt, sin, cos 
    111118 
    112119name = "elliptical_cylinder" 
    … …  
    125132              ["sld",         "1e-6/Ang^2", 4.0,   [-inf, inf], "sld",         "Cylinder scattering length density"], 
    126133              ["sld_solvent", "1e-6/Ang^2", 1.0,   [-inf, inf], "sld",         "Solvent scattering length density"], 
    127               ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "In plane angle"], 
    128               ["phi",         "degrees",    0,     [-360, 360], "orientation", "Out of plane angle"], 
    129               ["psi",         "degrees",    0,     [-360, 360], "orientation", "Major axis angle relative to Q"]] 
     134              ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "cylinder axis to beam angle"], 
     135              ["phi",         "degrees",    0,     [-360, 360], "orientation", "rotation about beam"], 
     136              ["psi",         "degrees",    0,     [-360, 360], "orientation", "rotation about cylinder axis"]] 
    130137 
    131138# pylint: enable=bad-whitespace, line-too-long 
    … …  
    149156                           + (length + radius) * (length + pi * radius)) 
    150157    return 0.5 * (ddd) ** (1. / 3.) 
     158q = 0.1 
     159# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
     160qx = q*cos(pi/6.0) 
     161qy = q*sin(pi/6.0) 
    151162 
    152163tests = [ 
    … …  
    158169      'sld_solvent':1.0, 'background':0.0}, 
    159170     0.001, 675.504402], 
     171#    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
    160172] 

  • sasmodels/models/fcc_paracrystal.c

    r4962519 r50beefe  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, cos_a1, cos_a2, cos_a3; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
     92    double q, zhat, yhat, xhat; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    9494 
    95     const double a1 = cos_a2 + cos_a3; 
    96     const double a2 = cos_a3 + cos_a1; 
    97     const double a3 = cos_a2 + cos_a1; 
     95    const double a1 = yhat + xhat; 
     96    const double a2 = xhat + zhat; 
     97    const double a3 = yhat + zhat; 
    9898    const double qd = 0.5*q*dnn; 
    9999    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 

  • sasmodels/models/fcc_paracrystal.py

    r925ad6e r69e1afc  
    76762D model computation. 
    7777 
    78 .. figure:: img/bcc_angle_definition.png 
     78.. figure:: img/parallelepiped_angle_definition.png 
    7979 
    80     Orientation of the crystal with respect to the scattering plane. 
     80    Orientation of the crystal with respect to the scattering plane, when  
     81    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8182 
    8283References 
    … …  
    9091""" 
    9192 
    92 from numpy import inf 
     93from numpy import inf, pi 
    9394 
    9495name = "fcc_paracrystal" 
    … …  
    110111              ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], 
    111112              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"], 
    112               ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 
    113               ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 
    114               ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 
     113              ["theta",       "degrees",    60,    [-360, 360], "orientation", "c axis to beam angle"], 
     114              ["phi",         "degrees",    60,    [-360, 360], "orientation", "rotation about beam"], 
     115              ["psi",         "degrees",    60,    [-360, 360], "orientation", "rotation about c axis"] 
    115116             ] 
    116117# pylint: enable=bad-whitespace, line-too-long 
    … …  
    128129            psi_pd=15, psi_pd_n=0, 
    129130           ) 
     131# april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     132q =4.*pi/220. 
     133tests = [ 
     134    [{ }, 
     135     [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]], 
     136     [{}, (-0.047,-0.007), 238.103096286], 
     137     [{}, (0.053,0.063), 0.863609587796 ], 
     138] 

  • sasmodels/models/hayter_msa.py

    ra807206 r5702f52  
    7474    ["radius_effective", "Ang", 20.75,   [0, inf],    "volume", "effective radius of charged sphere"], 
    7575    ["volfraction",   "None",     0.0192, [0, 0.74],   "", "volume fraction of spheres"], 
    76     ["charge",        "e",   19.0,    [0, inf],    "", "charge on sphere (in electrons)"], 
    77     ["temperature",   "K",  318.16,   [0, inf],    "", "temperature, in Kelvin, for Debye length calculation"], 
    78     ["concentration_salt",      "M",    0.0,    [-inf, inf], "", "conc of salt, moles/litre, 1:1 electolyte, for Debye length"], 
     76    ["charge",        "e",   19.0,    [0, 200],    "", "charge on sphere (in electrons)"], 
     77    ["temperature",   "K",  318.16,   [0, 450],    "", "temperature, in Kelvin, for Debye length calculation"], 
     78    ["concentration_salt",      "M",    0.0,    [0, inf], "", "conc of salt, moles/litre, 1:1 electolyte, for Debye length"], 
    7979    ["dielectconst",  "None",    71.08,   [-inf, inf], "", "dielectric constant (relative permittivity) of solvent, default water, for Debye length"] 
    8080    ] 

  • sasmodels/models/hollow_cylinder.py

    raea2e2a r9b79f29  
    6060""" 
    6161 
    62 from numpy import pi, inf 
     62from numpy import pi, inf, sin, cos 
    6363 
    6464name = "hollow_cylinder" 
    … …  
    8282    ["sld",         "1/Ang^2",  6.3, [-inf, inf], "sld",         "Cylinder sld"], 
    8383    ["sld_solvent", "1/Ang^2",  1,   [-inf, inf], "sld",         "Solvent sld"], 
    84     ["theta",       "degrees", 90,   [-360, 360], "orientation", "Theta angle"], 
    85     ["phi",         "degrees",  0,   [-360, 360], "orientation", "Phi angle"], 
     84    ["theta",       "degrees", 90,   [-360, 360], "orientation", "Cylinder axis to beam angle"], 
     85    ["phi",         "degrees",  0,   [-360, 360], "orientation", "Rotation about beam"], 
    8686    ] 
    8787# pylint: enable=bad-whitespace, line-too-long 
    … …  
    129129            theta_pd=10, theta_pd_n=5, 
    130130           ) 
    131  
     131q = 0.1 
     132# april 6 2017, rkh added a 2d unit test, assume correct! 
     133qx = q*cos(pi/6.0) 
     134qy = q*sin(pi/6.0) 
    132135# Parameters for unit tests 
    133136tests = [ 
    134137    [{}, 0.00005, 1764.926], 
    135138    [{}, 'VR', 1.8], 
    136     [{}, 0.001, 1756.76] 
    137     ] 
     139    [{}, 0.001, 1756.76], 
     140    [{}, (qx, qy), 2.36885476192  ], 
     141        ] 
     142del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/lib/sas_3j1x_x.c

    r473a9f1 reb2946f  
    4646double sas_3j1x_x(double q) 
    4747{ 
    48     if (q < SPH_J1C_CUTOFF) { 
     48    // 2017-05-18 PAK - support negative q 
     49    if (fabs(q) < SPH_J1C_CUTOFF) { 
    4950        const double q2 = q*q; 
    5051        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); 

  • sasmodels/models/lib/sas_J0.c

    rc8902ac reb2946f  
    236236        xx = x; 
    237237 
    238     if( x <= 2.0 ) { 
     238    // 2017-05-18 PAK - support negative x 
     239    if( xx <= 2.0 ) { 
    239240        z = xx * xx; 
    240         if( x < 1.0e-3 ) 
     241        if( xx < 1.0e-3 ) 
    241242            return( 1.0 - 0.25*z ); 
    242243 
    … …  
    245246    } 
    246247 
    247     q = 1.0/x; 
     248    q = 1.0/xx; 
    248249    w = sqrt(q); 
    249250 

  • sasmodels/models/lib/sas_J1.c

    r473a9f1 reb2946f  
    109109{ 
    110110 
    111     double w, z, p, q, xn; 
     111    double w, z, p, q, xn, abs_x, sign_x; 
    112112 
    113113    const double Z1 = 1.46819706421238932572E1; 
    … …  
    116116    const double SQ2OPI = 0.79788456080286535588; 
    117117 
    118     w = x; 
    119     if( x < 0 ) 
    120         w = -x; 
    121  
    122     if( w <= 5.0 ) { 
    123         z = x * x; 
     118    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x) 
     119    if (x < 0) { 
     120        abs_x = -x; 
     121        sign_x = -1.0; 
     122    } else { 
     123        abs_x = x; 
     124        sign_x = 1.0; 
     125    } 
     126 
     127    if( abs_x <= 5.0 ) { 
     128        z = abs_x * abs_x; 
    124129        w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
    125         w = w * x * (z - Z1) * (z - Z2); 
    126         return( w ); 
     130        w = w * abs_x * (z - Z1) * (z - Z2); 
     131        return( sign_x * w ); 
    127132    } 
    128133 
    129     w = 5.0/x; 
     134    w = 5.0/abs_x; 
    130135    z = w * w; 
    131  
    132136    p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 
    133137    q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 
    134  
    135     xn = x - THPIO4; 
     138    xn = abs_x - THPIO4; 
    136139 
    137140    double sn, cn; 
    … …  
    139142    p = p * cn - w * q * sn; 
    140143 
    141     return( p * SQ2OPI / sqrt(x) ); 
     144    return( sign_x * p * SQ2OPI / sqrt(abs_x) ); 
    142145} 
    143146 
    … …  
    179182    }; 
    180183 
    181 float cephes_j1f(float x) 
     184float cephes_j1f(float xx) 
    182185{ 
    183186 
    184     float xx, w, z, p, q, xn; 
     187    float x, w, z, p, q, xn; 
    185188 
    186189    const float Z1 = 1.46819706421238932572E1; 
    … …  
    188191 
    189192 
    190     xx = x; 
    191     if( xx < 0 ) 
    192         xx = -x; 
    193  
    194     if( xx <= 2.0 ) { 
    195         z = xx * xx; 
    196         p = (z-Z1) * xx * polevl( z, JPJ1, 4 ); 
    197         return( p ); 
     193    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x) 
     194    x = xx; 
     195    if( x < 0 ) 
     196        x = -xx; 
     197 
     198    if( x <= 2.0 ) { 
     199        z = x * x; 
     200        p = (z-Z1) * x * polevl( z, JPJ1, 4 ); 
     201        return( xx < 0. ? -p : p ); 
    198202    } 
    199203 
    … …  
    204208    w = q*q; 
    205209    xn = q * polevl( w, PH1J1, 7) - THPIO4F; 
    206     p = p * cos(xn + xx); 
    207  
    208     return(p); 
     210    p = p * cos(xn + x); 
     211 
     212    return( xx < 0. ? -p : p ); 
    209213} 
    210214#endif 

  • sasmodels/models/lib/sas_erf.c

    rb3796fa reb2946f  
    165165        // the erf function instead and z < 1.0. 
    166166        //return (1.0 - cephes_erf(a)); 
    167         z = x * x; 
    168         y = x * polevl(z, TD, 4) / p1evl(z, UD, 5); 
    169  
    170         return y; 
     167        // 2017-05-18 PAK - use erf(a) rather than erf(|a|) 
     168        z = a * a; 
     169        y = a * polevl(z, TD, 4) / p1evl(z, UD, 5); 
     170 
     171        return 1.0 - y; 
    171172    } 
    172173 
    … …  
    279280        //is explicit here for the case < 1.0 
    280281        //return (1.0 - sas_erf(a)); 
    281         z = x * x; 
    282         y = x * polevl( z, TF, 6 ); 
    283  
    284         return y; 
     282        // 2017-05-18 PAK - use erf(a) rather than erf(|a|) 
     283        z = a * a; 
     284        y = a * polevl( z, TF, 6 ); 
     285 
     286        return 1.0 - y; 
    285287    } 
    286288 

  • sasmodels/models/onion.py

    r925ad6e rbccb40f  
    11r""" 
    22This model provides the form factor, $P(q)$, for a multi-shell sphere where 
    3 the scattering length density (SLD) of the each shell is described by an 
     3the scattering length density (SLD) of each shell is described by an 
    44exponential, linear, or constant function. The form factor is normalized by 
    55the volume of the sphere where the SLD is not identical to the SLD of the 
    … …  
    142142 
    143143For $A = 0$, the exponential function has no dependence on the radius (so that 
    144 $\rho_\text{out}$ is ignored this case) and becomes flat. We set the constant 
     144$\rho_\text{out}$ is ignored in this case) and becomes flat. We set the constant 
    145145to $\rho_\text{in}$ for convenience, and thus the form factor contributed by 
    146146the shells is 
    … …  
    346346            # flat shell 
    347347            z.append(z_current + thickness[k]) 
    348             rho.append(sld_out[k]) 
     348            rho.append(sld_in[k]) 
    349349        else: 
    350350            # exponential shell 
    … …  
    357357                z.append(z_current+z_shell) 
    358358                rho.append(slope*exp(A[k]*z_shell/thickness[k]) + const) 
    359  
     359     
    360360    # add in the solvent 
    361361    z.append(z[-1]) 

  • sasmodels/models/parallelepiped.c

    r1e7b0db0 rd605080  
    6767    double psi) 
    6868{ 
    69     double q, cos_val_a, cos_val_b, cos_val_c; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
     69    double q, xhat, yhat, zhat; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    7171 
    72     const double siA = sas_sinx_x(0.5*q*length_a*cos_val_a); 
    73     const double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 
    74     const double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 
     72    const double siA = sas_sinx_x(0.5*length_a*q*xhat); 
     73    const double siB = sas_sinx_x(0.5*length_b*q*yhat); 
     74    const double siC = sas_sinx_x(0.5*length_c*q*zhat); 
    7575    const double V = form_volume(length_a, length_b, length_c); 
    7676    const double drho = (sld - solvent_sld); 

  • sasmodels/models/parallelepiped.py

    r8e68ea0 r34a9e4e  
    99---------- 
    1010 
    11 | 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`. 
     11 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`. 
    1414 
    1515.. _parallelepiped-image: 
    1616 
     17 
    1718.. figure:: img/parallelepiped_geometry.jpg 
    1819 
    … …  
    2122.. note:: 
    2223 
    23    The edge of the solid must satisfy the condition that $A < B < C$. 
    24    This requirement is not enforced in the model, so it is up to the 
    25    user to check this during the analysis. 
     24The three dimensions of the parallelepiped (strictly here a cuboid) may be given in  
     25$any$ size order. To avoid multiple fit solutions, especially 
     26with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may  
     27be a number of closely similar "best fits", so some trial and error, or fixing of some  
     28dimensions at expected values, may help. 
    2629 
    2730The 1D scattering intensity $I(q)$ is calculated as: 
    … …  
    6770    \mu &= qB 
    6871 
    69  
    7072The scattering intensity per unit volume is returned in units of |cm^-1|. 
    7173 
    7274NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    73 the averaged effective radius $(=\sqrt{A B / \pi})$ and 
     75the averaged effective radius, after appropriately sorting the three 
     76dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
    7477length $(= C)$ values, and used as the effective radius for 
    7578$S(q)$ when $P(q) \cdot S(q)$ is applied. 
    … …  
    102105.. _parallelepiped-orientation: 
    103106 
    104 .. figure:: img/parallelepiped_angle_definition.jpg 
    105  
    106     Definition of the angles for oriented parallelepipeds. 
    107  
    108 .. figure:: img/parallelepiped_angle_projection.jpg 
    109  
    110     Examples of the angles for oriented parallelepipeds against the 
     107.. figure:: img/parallelepiped_angle_definition.png 
     108 
     109    Definition of the angles for oriented parallelepiped, shown with $A<B<C$. 
     110 
     111.. figure:: img/parallelepiped_angle_projection.png 
     112 
     113    Examples of the angles for an oriented parallelepiped against the 
    111114    detector plane. 
    112115 
     116On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
     117appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces.  
     118(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is  
     119about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to  
     120understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation  
     121distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)  
     122 
     123     
    113124For a given orientation of the parallelepiped, the 2D form factor is 
    114125calculated as 
    … …  
    116127.. math:: 
    117128 
    118     P(q_x, q_y) = \left[\frac{\sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2 
    119                   \left[\frac{\sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2 
    120                   \left[\frac{\sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2 
     129    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 
     130                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 
     131                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 
    121132 
    122133with 
    … …  
    154165angles. 
    155166 
    156 This model is based on form factor calculations implemented in a c-library 
    157 provided by the NIST Center for Neutron Research (Kline, 2006). 
    158167 
    159168References 
    … …  
    163172 
    164173R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     174 
     175Authorship and Verification 
     176---------------------------- 
     177 
     178* **Author:** This model is based on form factor calculations implemented 
     179    in a c-library provided by the NIST Center for Neutron Research (Kline, 2006). 
     180* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017 
     181* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017 
     182 
    165183""" 
    166184 
    167185import numpy as np 
    168 from numpy import pi, inf, sqrt 
     186from numpy import pi, inf, sqrt, sin, cos 
    169187 
    170188name = "parallelepiped" 
    … …  
    180198            mu = q*B 
    181199        V: Volume of the rectangular parallelepiped 
    182         alpha: angle between the long axis of the  
     200        alpha: angle between the long axis of the 
    183201            parallelepiped and the q-vector for 1D 
    184202""" 
    … …  
    196214              ["length_c", "Ang", 400, [0, inf], "volume", 
    197215               "Larger side of the parallelepiped"], 
    198               ["theta", "degrees", 60, [-inf, inf], "orientation", 
    199                "In plane angle"], 
    200               ["phi", "degrees", 60, [-inf, inf], "orientation", 
    201                "Out of plane angle"], 
    202               ["psi", "degrees", 60, [-inf, inf], "orientation", 
    203                "Rotation angle around its own c axis against q plane"], 
     216              ["theta", "degrees", 60, [-360, 360], "orientation", 
     217               "c axis to beam angle"], 
     218              ["phi", "degrees", 60, [-360, 360], "orientation", 
     219               "rotation about beam"], 
     220              ["psi", "degrees", 60, [-360, 360], "orientation", 
     221               "rotation about c axis"], 
    204222             ] 
    205223 
    … …  
    208226def ER(length_a, length_b, length_c): 
    209227    """ 
    210         Return effective radius (ER) for P(q)*S(q) 
     228    Return effective radius (ER) for P(q)*S(q) 
    211229    """ 
    212  
     230    # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller 
     231    # or much larger 
     232    abc = np.vstack((length_a, length_b, length_c)) 
     233    abc = np.sort(abc, axis=0) 
     234    selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 
     235    length = np.where(selector, abc[0], abc[2]) 
    213236    # surface average radius (rough approximation) 
    214     surf_rad = sqrt(length_a * length_b / pi) 
    215  
    216     ddd = 0.75 * surf_rad * (2 * surf_rad * length_c + (length_c + surf_rad) * (length_c + pi * surf_rad)) 
     237    radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 
     238 
     239    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 
    217240    return 0.5 * (ddd) ** (1. / 3.) 
    218241 
    … …  
    230253            phi_pd=10, phi_pd_n=1, 
    231254            psi_pd=10, psi_pd_n=10) 
    232  
    233 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     255# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     256qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    234257tests = [[{}, 0.2, 0.17758004974], 
    235258         [{}, [0.2], [0.17758004974]], 
    236          [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014], 
    237          [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]], 
     259         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475], 
     260         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]], 
    238261        ] 
    239262del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/sc_paracrystal.c

    r4962519 r50beefe  
    111111          double psi) 
    112112{ 
    113     double q, cos_a1, cos_a2, cos_a3; 
    114     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
     113    double q, zhat, yhat, xhat; 
     114    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    115115 
    116116    const double qd = q*dnn; 
    … …  
    118118    const double tanh_qd = tanh(arg); 
    119119    const double cosh_qd = cosh(arg); 
    120     const double Zq = tanh_qd/(1. - cos(qd*cos_a1)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd); 
     120    const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 
     121                    * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
     122                    * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
    123123 
    124124    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 

  • sasmodels/models/sc_paracrystal.py

    r925ad6e r69e1afc  
    8383model computation. 
    8484 
    85 .. figure:: img/sc_crystal_angle_definition.jpg 
     85.. figure:: img/parallelepiped_angle_definition.png 
     86 
     87    Orientation of the crystal with respect to the scattering plane, when 
     88    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8689 
    8790Reference 
    … …  
    127130              ["sld",  "1e-6/Ang^2",         3.0, [0.0, inf],  "sld",         "Sphere scattering length density"], 
    128131              ["sld_solvent", "1e-6/Ang^2",  6.3, [0.0, inf],  "sld",         "Solvent scattering length density"], 
    129               ["theta",       "degrees",     0.0, [-inf, inf], "orientation", "Orientation of the a1 axis w/respect incoming beam"], 
    130               ["phi",         "degrees",     0.0, [-inf, inf], "orientation", "Orientation of the a2 in the plane of the detector"], 
    131               ["psi",         "degrees",     0.0, [-inf, inf], "orientation", "Orientation of the a3 in the plane of the detector"], 
     132              ["theta",       "degrees",    0,    [-360, 360], "orientation", "c axis to beam angle"], 
     133              ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
     134              ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about c axis"] 
    132135             ] 
    133136# pylint: enable=bad-whitespace, line-too-long 
    … …  
    146149 
    147150tests = [ 
    148     # Accuracy tests based on content in test/utest_extra_models.py 
     151    # Accuracy tests based on content in test/utest_extra_models.py, 2d tests added April 10, 2017 
    149152    [{}, 0.001, 10.3048], 
    150153    [{}, 0.215268, 0.00814889], 
    151     [{}, (0.414467), 0.001313289] 
     154    [{}, (0.414467), 0.001313289], 
     155    [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ], 
     156    [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ] 
    152157    ] 
    153158 

  • sasmodels/models/spherical_sld.py

    r925ad6e r63a7fe8  
    199199category = "shape:sphere" 
    200200 
    201 SHAPES = [["erf(|nu|*z)", "Rpow(z^|nu|)", "Lpow(z^|nu|)", 
    202            "Rexp(-|nu|z)", "Lexp(-|nu|z)"]] 
     201SHAPES = ["erf(|nu|*z)", "Rpow(z^|nu|)", "Lpow(z^|nu|)", 
     202          "Rexp(-|nu|z)", "Lexp(-|nu|z)"] 
    203203 
    204204# pylint: disable=bad-whitespace, line-too-long 
    … …  
    209209              ["thickness[n_shells]",  "Ang",        100.0,  [0, inf],       "volume", "thickness shell"], 
    210210              ["interface[n_shells]",  "Ang",        50.0,   [0, inf],       "volume", "thickness of the interface"], 
    211               ["shape[n_shells]",      "",           0,      SHAPES,         "", "interface shape"], 
     211              ["shape[n_shells]",      "",           0,      [SHAPES],       "", "interface shape"], 
    212212              ["nu[n_shells]",         "",           2.5,    [0, inf],       "", "interface shape exponent"], 
    213213              ["n_steps",              "",           35,     [0, inf],       "", "number of steps in each interface (must be an odd integer)"], 

  • sasmodels/models/stacked_disks.py

    rb7e8b94 r1b85b55  
    5858and $\sigma_d$ = the Gaussian standard deviation of the d-spacing (*sigma_d*). 
    5959Note that $D\cos(\alpha)$ is the component of $D$ parallel to $q$ and the last 
    60 term in the equation above is effectively a Debye-Waller factor term.  
     60term in the equation above is effectively a Debye-Waller factor term. 
    6161 
    6262.. note:: 
    … …  
    7777the axis of the cylinder using two angles $\theta$ and $\varphi$. 
    7878 
    79 .. figure:: img/cylinder_angle_definition.jpg 
     79.. figure:: img/cylinder_angle_definition.png 
    8080 
    8181    Examples of the angles against the detector plane. 
    … …  
    103103""" 
    104104 
    105 from numpy import inf 
     105from numpy import inf, sin, cos, pi 
    106106 
    107107name = "stacked_disks" 
    … …  
    131131    ["sld_layer",   "1e-6/Ang^2",  0.0, [-inf, inf], "sld",         "Layer scattering length density"], 
    132132    ["sld_solvent", "1e-6/Ang^2",  5.0, [-inf, inf], "sld",         "Solvent scattering length density"], 
    133     ["theta",       "degrees",     0,   [-inf, inf], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"], 
    134     ["phi",         "degrees",     0,   [-inf, inf], "orientation", "Orientation of the stacked disk in the plane of the detector"], 
     133    ["theta",       "degrees",     0,   [-360, 360], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"], 
     134    ["phi",         "degrees",     0,   [-360, 360], "orientation", "Rotation about beam"], 
    135135    ] 
    136136# pylint: enable=bad-whitespace, line-too-long 
    … …  
    152152# After redefinition of spherical coordinates - 
    153153# tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
    154 # but should not matter here as so far all the tests are 1D not 2D 
     154q = 0.1 
     155# april 6 2017, rkh added a 2d unit test, assume correct! 
     156qx = q*cos(pi/6.0) 
     157qy = q*sin(pi/6.0) 
     158# Accuracy tests based on content in test/utest_extra_models.py. 
     159# Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; 
     160# for which alas q=0.001 values seem closer to n_stacked=1 not 5, 
     161# changed assuming my 4.1 code OK, RKH 
    155162tests = [ 
    156 # Accuracy tests based on content in test/utest_extra_models.py. 
    157 # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; for which alas q=0.001 values seem closer to n_stacked=1 not 5, changed assuming my 4.1 code OK, RKH 
    158     [{'thick_core': 10.0, 
    159       'thick_layer': 15.0, 
    160       'radius': 3000.0, 
    161       'n_stacking': 1.0, 
    162       'sigma_d': 0.0, 
    163       'sld_core': 4.0, 
    164       'sld_layer': -0.4, 
    165       'solvent_sd': 5.0, 
     163    [{'thick_core': 10.0, 
     164      'thick_layer': 15.0, 
     165      'radius': 3000.0, 
     166      'n_stacking': 1.0, 
     167      'sigma_d': 0.0, 
     168      'sld_core': 4.0, 
     169      'sld_layer': -0.4, 
     170      'sld_solvent': 5.0, 
    166171      'theta': 90.0, 
    167172      'phi': 0.0, 
    … …  
    176181      'sld_core': 4.0, 
    177182      'sld_layer': -0.4, 
    178       'solvent_sd': 5.0, 
     183      'sld_solvent': 5.0, 
    179184      'theta': 90.0, 
    180185      'phi': 0.0, 
    … …  
    186191    [{'thick_core': 10.0, 
    187192      'thick_layer': 15.0, 
    188       'radius': 3000.0, 
     193      'radius': 100.0, 
    189194      'n_stacking': 5, 
    190195      'sigma_d': 0.0, 
    191196      'sld_core': 4.0, 
    192197      'sld_layer': -0.4, 
    193       'solvent_sd': 5.0, 
     198      'sld_solvent': 5.0, 
     199      'theta': 90.0, 
     200      'phi': 20.0, 
     201      'scale': 0.01, 
     202      'background': 0.001, 
     203     }, (qx, qy), 0.0491167089952], 
     204    [{'thick_core': 10.0, 
     205      'thick_layer': 15.0, 
     206      'radius': 3000.0, 
     207      'n_stacking': 5, 
     208      'sigma_d': 0.0, 
     209      'sld_core': 4.0, 
     210      'sld_layer': -0.4, 
     211      'sld_solvent': 5.0, 
    194212      'theta': 90.0, 
    195213      'phi': 0.0, 
    … …  
    207225      'sld_core': 4.0, 
    208226      'sld_layer': -0.4, 
    209       'solvent_sd': 5.0, 
     227      'sld_solvent': 5.0, 
    210228      'theta': 90.0, 
    211229      'phi': 0.0, 
    … …  
    222240      'sld_core': 4.0, 
    223241      'sld_layer': -0.4, 
    224       'solvent_sd': 5.0, 
     242      'sld_solvent': 5.0, 
    225243      'theta': 90.0, 
    226244      'phi': 0.0, 
    … …  
    228246      'background': 0.001, 
    229247     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
    230  
    231     [{'thick_core': 10.0, 
    232       'thick_layer': 15.0, 
    233       'radius': 3000.0, 
    234       'n_stacking': 1.0, 
    235       'sigma_d': 0.0, 
    236       'sld_core': 4.0, 
    237       'sld_layer': -0.4, 
    238      'solvent_sd': 5.0, 
     248    [{'thick_core': 10.0, 
     249      'thick_layer': 15.0, 
     250      'radius': 3000.0, 
     251      'n_stacking': 1.0, 
     252      'sigma_d': 0.0, 
     253      'sld_core': 4.0, 
     254      'sld_layer': -0.4, 
     255      'sld_solvent': 5.0, 
     256      'theta': 90.0, 
     257      'phi': 20.0, 
     258      'scale': 0.01, 
     259      'background': 0.001, 
     260# 2017-05-18 PAK temporarily suppress output of qx,qy test; j1 is not accurate for large qr 
     261#     }, (qx, qy), 0.0341738733124], 
     262     }, (qx, qy), None], 
     263 
     264    [{'thick_core': 10.0, 
     265      'thick_layer': 15.0, 
     266      'radius': 3000.0, 
     267      'n_stacking': 1.0, 
     268      'sigma_d': 0.0, 
     269      'sld_core': 4.0, 
     270      'sld_layer': -0.4, 
     271      'sld_solvent': 5.0, 
    239272      'theta': 90.0, 
    240273      'phi': 0.0, 
    … …  
    243276     }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
    244277    ] 
    245 # 11Jan2017   RKH checking unit test agai, note they are all 1D, no 2D 
    246  
     278# 11Jan2017   RKH checking unit test again, note they are all 1D, no 2D 

  • sasmodels/models/triaxial_ellipsoid.c

    r925ad6e r68dd6a9  
    2020    double radius_polar) 
    2121{ 
    22     double sn, cn; 
    23     // translate a point in [-1,1] to a point in [0, 1] 
    24     const double zm = 0.5; 
    25     const double zb = 0.5; 
     22    const double pa = square(radius_equat_minor/radius_equat_major) - 1.0; 
     23    const double pc = square(radius_polar/radius_equat_major) - 1.0; 
     24    // translate a point in [-1,1] to a point in [0, pi/2] 
     25    const double zm = M_PI_4; 
     26    const double zb = M_PI_4; 
    2627    double outer = 0.0; 
    2728    for (int i=0;i<76;i++) { 
    28         //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    29         const double x = 0.5*(Gauss76Z[i] + 1.0); 
    30         SINCOS(M_PI_2*x, sn, cn); 
    31         const double acosx2 = radius_equat_minor*radius_equat_minor*cn*cn; 
    32         const double bsinx2 = radius_equat_major*radius_equat_major*sn*sn; 
    33         const double c2 = radius_polar*radius_polar; 
     29        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 
     30        const double phi = Gauss76Z[i]*zm + zb; 
     31        const double pa_sinsq_phi = pa*square(sin(phi)); 
    3432 
    3533        double inner = 0.0; 
     34        const double um = 0.5; 
     35        const double ub = 0.5; 
    3636        for (int j=0;j<76;j++) { 
    37             const double ysq = square(Gauss76Z[j]*zm + zb); 
    38             const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq); 
    39             const double fq = sas_3j1x_x(t); 
    40             inner += Gauss76Wt[j] * fq * fq ; 
     37            // translate a point in [-1,1] to a point in [0, 1] 
     38            const double usq = square(Gauss76Z[j]*um + ub); 
     39            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
     40            const double fq = sas_3j1x_x(q*r); 
     41            inner += Gauss76Wt[j] * fq * fq; 
    4142        } 
    42         outer += Gauss76Wt[i] * 0.5 * inner; 
     43        outer += Gauss76Wt[i] * inner;  // correcting for dx later 
    4344    } 
    44     // translate dx in [-1,1] to dx in [lower,upper] 
    45     const double fqsq = outer*zm; 
     45    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
     46    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    4647    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    4748    return 1.0e-4 * s * s * fqsq; 
    … …  
    5859    double psi) 
    5960{ 
    60     double q, calpha, cmu, cnu; 
    61     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 
     61    double q, xhat, yhat, zhat; 
     62    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    6263 
    63     const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu 
    64                           + radius_equat_major*radius_equat_major*cmu*cmu 
    65                           + radius_polar*radius_polar*calpha*calpha); 
    66     const double fq = sas_3j1x_x(t); 
     64    const double r = sqrt(square(radius_equat_minor*xhat) 
     65                          + square(radius_equat_major*yhat) 
     66                          + square(radius_polar*zhat)); 
     67    const double fq = sas_3j1x_x(q*r); 
    6768    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    6869 

  • sasmodels/models/triaxial_ellipsoid.py

    r925ad6e r34a9e4e  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 All three axes are of different lengths with $R_a \leq R_b \leq R_c$ 
    5 **Users should maintain this inequality for all calculations**. 
     4Definition 
     5---------- 
     6 
     7.. figure:: img/triaxial_ellipsoid_geometry.jpg 
     8 
     9    Ellipsoid with $R_a$ as *radius_equat_minor*, $R_b$ as *radius_equat_major* 
     10    and $R_c$ as *radius_polar*. 
     11 
     12Given an ellipsoid 
    613 
    714.. math:: 
    815 
    9     P(q) = \text{scale} V \left< F^2(q) \right> + \text{background} 
     16    \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1 
    1017 
    11 where the volume $V = 4/3 \pi R_a R_b R_c$, and the averaging 
    12 $\left<\ldots\right>$ is applied over all orientations for 1D. 
    13  
    14 .. figure:: img/triaxial_ellipsoid_geometry.jpg 
    15  
    16     Ellipsoid schematic. 
    17  
    18 Definition 
    19 ---------- 
    20  
    21 The form factor calculated is 
     18the scattering for randomly oriented particles is defined by the average over 
     19all orientations $\Omega$ of: 
    2220 
    2321.. math:: 
    2422 
    25     P(q) = \frac{\text{scale}}{V}\int_0^1\int_0^1 
    26         \Phi^2(qR_a^2\cos^2( \pi x/2) + qR_b^2\sin^2(\pi y/2)(1-y^2) + R_c^2y^2) 
    27         dx dy 
     23    P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega\Phi^2(qr)\,d\Omega 
     24           + \text{background} 
    2825 
    2926where 
    … …  
    3128.. math:: 
    3229 
    33     \Phi(u) = 3 u^{-3} (\sin u - u \cos u) 
     30    \Phi(qr) &= 3 j_1(qr)/qr = 3 (\sin qr - qr \cos qr)/(qr)^3 \\ 
     31    r^2 &= R_a^2e^2 + R_b^2f^2 + R_c^2g^2 \\ 
     32    V &= \tfrac{4}{3} \pi R_a R_b R_c 
    3433 
     34The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$, 
     35$Y$ and $Z$ respectively.  Keeping the orientation fixed at the canonical 
     36axes, we can integrate over the incident direction using polar angle 
     37$-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$ 
     38(as defined in ref [1]), 
     39 
     40 .. math:: 
     41 
     42     \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr) 
     43                                                \cos \gamma\,d\gamma d\phi 
     44 
     45with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$. 
     46A little algebra yields 
     47 
     48.. math:: 
     49 
     50    r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma) 
     51 
     52for 
     53 
     54.. math:: 
     55 
     56    p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1 
     57 
     58Due to symmetry, the ranges can be restricted to a single quadrant 
     59$0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting 
     60integral by 8. The computation is done using the substitution $u = \sin\gamma$, 
     61$du = \cos\gamma\,d\gamma$, giving 
     62 
     63.. math:: 
     64 
     65    \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\ 
     66    r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2) 
     67 
     68Though for convenience we describe the three radii of the ellipsoid as equatorial 
     69and polar, they may be given in $any$ size order. To avoid multiple solutions, especially 
     70with Monte-Carlo fit methods, it may be advisable to restrict their ranges. For typical 
     71small angle diffraction situations there may be a number of closely similar "best fits", 
     72so some trial and error, or fixing of some radii at expected values, may help. 
     73     
    3574To provide easy access to the orientation of the triaxial ellipsoid, 
    3675we define the axis of the cylinder using the angles $\theta$, $\phi$ 
    37 and $\psi$. These angles are defined on 
    38 :numref:`triaxial-ellipsoid-angles` . 
    39 The angle $\psi$ is the rotational angle around its own $c$ axis 
    40 against the $q$ plane. For example, $\psi = 0$ when the 
    41 $a$ axis is parallel to the $x$ axis of the detector. 
     76and $\psi$. These angles are defined analogously to the elliptical_cylinder below, note that 
     77angle $\phi$ is now NOT the same as in the equations above. 
     78 
     79.. figure:: img/elliptical_cylinder_angle_definition.png 
     80 
     81    Definition of angles for oriented triaxial ellipsoid, where radii are for illustration here  
     82    $a < b << c$ and angle $\Psi$ is a rotation around the axis of the particle. 
     83 
     84For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     85see the :ref:`elliptical-cylinder` model for further information. 
    4286 
    4387.. _triaxial-ellipsoid-angles: 
    4488 
    45 .. figure:: img/triaxial_ellipsoid_angle_projection.jpg 
     89.. figure:: img/triaxial_ellipsoid_angle_projection.png 
    4690 
    47     The angles for oriented ellipsoid. 
     91    Some examples for an oriented triaxial ellipsoid. 
    4892 
    4993The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$. 
    5094 
    51 The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the 
     95The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the 
    5296parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major 
    5397equatorial radius, and $R_c$ is the polar radius of the ellipsoid. 
    5498 
    5599NB: The 2nd virial coefficient of the triaxial solid ellipsoid is 
    56 calculated based on the polar radius $R_p = R_c$ and equatorial 
    57 radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for 
     100calculated after sorting the three radii to give the most appropriate 
     101prolate or oblate form, from the new polar radius $R_p = R_c$ and effective equatorial 
     102radius,  $R_e = \sqrt{R_a R_b}$, to then be used as the effective radius for 
    58103$S(q)$ when $P(q) \cdot S(q)$ is applied. 
    59104 
    … …  
    69114---------- 
    70115 
    71 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
    72 and Neutron Scattering*, Plenum, New York, 1987. 
     116[1] Finnigan, J.A., Jacobs, D.J., 1971. 
     117*Light scattering by ellipsoidal particles in solution*, 
     118J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 
     119 
     120Authorship and Verification 
     121---------------------------- 
     122 
     123* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     124* **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017 
     125* **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017 
    73126""" 
    74127 
    75 from numpy import inf 
     128from numpy import inf, sin, cos, pi 
    76129 
    77130name = "triaxial_ellipsoid" 
    78131title = "Ellipsoid of uniform scattering length density with three independent axes." 
    79132 
    80 description = """\ 
    81 Note: During fitting ensure that the inequality ra<rb<rc is not 
    82         violated. Otherwise the calculation will 
    83         not be correct. 
     133description = """ 
     134   Triaxial ellipsoid - see main documentation. 
    84135""" 
    85136category = "shape:ellipsoid" 
    … …  
    91142               "Solvent scattering length density"], 
    92143              ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 
    93                "Minor equatorial radius"], 
     144               "Minor equatorial radius, Ra"], 
    94145              ["radius_equat_major", "Ang", 400, [0, inf], "volume", 
    95                "Major equatorial radius"], 
     146               "Major equatorial radius, Rb"], 
    96147              ["radius_polar", "Ang", 10, [0, inf], "volume", 
    97                "Polar radius"], 
    98               ["theta", "degrees", 60, [-inf, inf], "orientation", 
    99                "In plane angle"], 
    100               ["phi", "degrees", 60, [-inf, inf], "orientation", 
    101                "Out of plane angle"], 
    102               ["psi", "degrees", 60, [-inf, inf], "orientation", 
    103                "Out of plane angle"], 
     148               "Polar radius, Rc"], 
     149              ["theta", "degrees", 60, [-360, 360], "orientation", 
     150               "polar axis to beam angle"], 
     151              ["phi", "degrees", 60, [-360, 360], "orientation", 
     152               "rotation about beam"], 
     153              ["psi", "degrees", 60, [-360, 360], "orientation", 
     154               "rotation about polar axis"], 
    104155             ] 
    105156 
    … …  
    108159def ER(radius_equat_minor, radius_equat_major, radius_polar): 
    109160    """ 
    110         Returns the effective radius used in the S*P calculation 
     161    Returns the effective radius used in the S*P calculation 
    111162    """ 
    112163    import numpy as np 
    113164    from .ellipsoid import ER as ellipsoid_ER 
    114     return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major)) 
     165 
     166    # now that radii can be in any size order, radii need sorting a,b,c 
     167    # where a~b and c is either much smaller or much larger 
     168    radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar)) 
     169    radii = np.sort(radii, axis=0) 
     170    selector = (radii[1] - radii[0]) > (radii[2] - radii[1]) 
     171    polar = np.where(selector, radii[0], radii[2]) 
     172    equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2])) 
     173    return ellipsoid_ER(polar, equatorial) 
    115174 
    116175demo = dict(scale=1, background=0, 
    … …  
    124183            phi_pd=15, phi_pd_n=1, 
    125184            psi_pd=15, psi_pd_n=1) 
     185 
     186q = 0.1 
     187# april 6 2017, rkh add unit tests 
     188#     NOT compared with any other calc method, assume correct! 
     189# check 2d test after pull #890 
     190qx = q*cos(pi/6.0) 
     191qy = q*sin(pi/6.0) 
     192tests = [[{}, 0.05, 24.8839548033], 
     193        [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266 ], 
     194        ] 
     195del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/two_power_law.py

    r40a87fa rbb4b509  
    120120     }, 0.150141, 0.125945], 
    121121 
    122     [{'coeffcent_1':    1.0, 
     122    [{'coefficent_1':    1.0, 
    123123      'crossover':  0.04, 
    124124      'power_1':    1.0, 
    … …  
    127127     }, 0.442528, 0.00166884], 
    128128 
    129     [{'coeffcent_1':    1.0, 
     129    [{'coefficent_1':    1.0, 
    130130      'crossover':  0.04, 
    131131      'power_1':    1.0, 

  • sasmodels/rst2html.py

    r0890871 rf2f5413  
    2929from docutils.writers.html4css1 import HTMLTranslator 
    3030from docutils.nodes import SkipNode 
    31  
    32 def wxview(html, url="", size=(850, 540)): 
    33     import wx 
    34     from wx.html2 import WebView 
    35     frame = wx.Frame(None, -1, size=size) 
    36     view = WebView.New(frame) 
    37     view.SetPage(html, url) 
    38     frame.Show() 
    39     return frame 
    40  
    41 def view_rst(filename): 
    42     from os.path import expanduser 
    43     with open(expanduser(filename)) as fid: 
    44         rst = fid.read() 
    45     html = rst2html(rst) 
    46     wxview(html) 
    4731 
    4832def rst2html(rst, part="whole", math_output="mathjax"): 
    … …  
    155139    assert replace_dollar(u"a (again $in parens$) a") == u"a (again :math:`in parens`) a" 
    156140 
    157 def view_rst_app(filename): 
     141def load_rst_as_html(filename): 
     142    from os.path import expanduser 
     143    with open(expanduser(filename)) as fid: 
     144        rst = fid.read() 
     145    html = rst2html(rst) 
     146    return html 
     147 
     148def wxview(html, url="", size=(850, 540)): 
     149    import wx 
     150    from wx.html2 import WebView 
     151    frame = wx.Frame(None, -1, size=size) 
     152    view = WebView.New(frame) 
     153    view.SetPage(html, url) 
     154    frame.Show() 
     155    return frame 
     156 
     157def view_html_wxapp(html, url=""): 
    158158    import wx  # type: ignore 
    159159    app = wx.App() 
    160     view_rst(filename) 
     160    frame = wxview(html, url) 
    161161    app.MainLoop() 
    162162 
     163def view_url_wxapp(url): 
     164    import wx  # type: ignore 
     165    from wx.html2 import WebView 
     166    app = wx.App() 
     167    frame = wx.Frame(None, -1, size=(850, 540)) 
     168    view = WebView.New(frame) 
     169    view.LoadURL(url) 
     170    frame.Show() 
     171    app.MainLoop() 
     172 
     173def qtview(html, url=""): 
     174    try: 
     175        from PyQt5.QtWebKitWidgets import QWebView 
     176        from PyQt5.QtCore import QUrl 
     177    except ImportError: 
     178        from PyQt4.QtWebkit import QWebView 
     179        from PyQt4.QtCore import QUrl 
     180    helpView = QWebView() 
     181    helpView.setHtml(html, QUrl(url)) 
     182    helpView.show() 
     183    return helpView 
     184 
     185def view_html_qtapp(html, url=""): 
     186    import sys 
     187    try: 
     188        from PyQt5.QtWidgets import QApplication 
     189    except ImportError: 
     190        from PyQt4.QtGui import QApplication 
     191    app = QApplication([]) 
     192    frame = qtview(html, url) 
     193    sys.exit(app.exec_()) 
     194 
     195def view_url_qtapp(url): 
     196    import sys 
     197    try: 
     198        from PyQt5.QtWidgets import QApplication 
     199    except ImportError: 
     200        from PyQt4.QtGui import QApplication 
     201    app = QApplication([]) 
     202    try: 
     203        from PyQt5.QtWebKitWidgets import QWebView 
     204        from PyQt5.QtCore import QUrl 
     205    except ImportError: 
     206        from PyQt4.QtWebkit import QWebView 
     207        from PyQt4.QtCore import QUrl 
     208    frame = QWebView() 
     209    frame.load(QUrl(url)) 
     210    frame.show() 
     211    sys.exit(app.exec_()) 
     212 
     213def view_help(filename, qt=False): 
     214    import os 
     215    url="file:///"+os.path.abspath(filename).replace("\\","/") 
     216    if filename.endswith('.rst'): 
     217        html = load_rst_as_html(filename) 
     218        if qt: 
     219            view_html_qtapp(html, url) 
     220        else: 
     221            view_html_wxapp(html, url) 
     222    else: 
     223        if qt: 
     224            view_url_qtapp(url) 
     225        else: 
     226            view_url_wxapp(url) 
    163227 
    164228if __name__ == "__main__": 
    165229    import sys 
    166     view_rst_app(sys.argv[1]) 
    167  
     230    view_help(sys.argv[1], qt=True) 
     231 
Note: See TracChangeset for help on using the changeset viewer.