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


Ignore:
Files:
15 added
14 deleted
50 edited

Legend:

Unmodified
Added
Removed

  • .travis.yml

    r24cd982 rbec0e75  
     1# Test Travis CL 
     2 
    13language: python 
    2  
    3 sudo:  false 
    4  
    5 matrix: 
    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  
     4python: 
     5  - "2.7" 
    266# whitelist 
    277branches: 
    288  only: 
    299    - master 
    30  
    31 addons: 
    32   apt: 
    33     packages: 
    34       opencl-headers 
    35  
     10# command to install dependencies 
     11virtualenv: 
     12  system_site_packages: true 
    3613before_install: 
    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; 
     14  - sudo apt-get update; 
     15  - sudo apt-get install python-numpy python-scipy 
     16#  - sudo apt-get install opencl-headers 
    5917 
    6018install: 
     19#  - pip install pyopencl 
    6120  - pip install bumps 
    6221  - pip install unittest-xml-reporting 
    63  
    6422script: 
    65 - export WORKSPACE=/home/travis/build/SasView/sasmodels/ 
    66 - python -m sasmodels.model_test dll all 
    67  
    68 notifications: 
    69   slack: 
    70     secure: xNAUeSu1/it/x9Q2CSg79aw1LLc7d6mLpcqSCTeKROp71RhkFf8VjJnJm/lEbKHNC8yj5H9UHrz5DmzwJzI+6oMt4NdEeS6WvGhwGY/wCt2IcJKxw0vj1DAU04qFMS041Khwclo6jIqm76DloinXvmvsS+K/nSyQkF7q4egSlwA= 
     23  - export WORKSPACE=/home/travis/build/SasView/sasmodels/ 
     24  - python -m sasmodels.model_test dll all 

  • LICENSE.txt

    rf68e2a5 rc8de1bd  
    1 Copyright (c) 2009-2017, SasView Developers 
    2 All rights reserved. 
     1This program is in the public domain. 
    32 
    4 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 
    5  
    6 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
    7  
    8 2. 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  
    10 3. 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  
    12 THIS 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. 
     3Individual files may be copyright different authors, with licences that 
     4allow modification and redistribution in source or binary form with or 
     5without modification, and a disclaimer of warranty. 

  • example/oriented_usans.py

    r1cd24b4 rea75043  
    1919    scale=0.08, background=0, 
    2020    sld=.291, sld_solvent=7.105, 
    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, 
     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, 
    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.radius_polar.range(1000, 10000) 
    29 model.radius_equatorial.range(1000, 10000) 
     28model.r_polar.range(1000, 10000) 
     29model.r_equatorial.range(1000, 10000) 
    3030model.theta.range(0, 360) 
    3131model.phi.range(0, 360) 

  • explore/angular_pd.py

    r8267e0b r12eb36b  
    4747 
    4848def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    49     theta_center = radians(90-theta) 
     49    theta_center = radians(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=90-theta if use_new else theta, azim=phi) 
     139            ax.view_init(elev=theta, azim=phi) 
    140140        plt.gcf().canvas.draw() 
    141141 

  • sasmodels/compare.py

    r630156b rf72d70a  
    3030 
    3131import sys 
    32 import os 
    3332import math 
    3433import datetime 
    … …  
    4140from . import kerneldll 
    4241from . import exception 
    43 from .data import plot_theory, empty_data1D, empty_data2D, load_data 
     42from .data import plot_theory, empty_data1D, empty_data2D 
    4443from .direct_model import DirectModel 
    4544from .convert import revert_name, revert_pars, constrain_new_to_old 
    … …  
    7372    -1d*/-2d computes 1d or 2d data 
    7473    -preset*/-random[=seed] preset or random parameters 
    75     -mono*/-poly force monodisperse or allow polydisperse demo parameters 
     74    -mono/-poly* force monodisperse/polydisperse 
    7675    -magnetic/-nonmagnetic* suppress magnetism 
    7776    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    … …  
    8483    -edit starts the parameter explorer 
    8584    -default/-demo* use demo vs default parameters 
    86     -help/-html shows the model docs instead of running the model 
     85    -html shows the model docs instead of running the model 
    8786    -title="note" adds note to the plot title, after the model name 
    88     -data="path" uses q, dq from the data file 
    8987 
    9088Any two calculation engines can be selected for comparison: 
    … …  
    546544    'f32': '32', 
    547545    'f64': '64', 
    548     'float16': '16', 
    549     'float32': '32', 
    550     'float64': '64', 
    551     'float128': '128', 
    552546    'longdouble': '128', 
    553547} 
    … …  
    753747    comp = opts['engines'][1] if have_comp else None 
    754748    data = opts['data'] 
    755     use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    756749 
    757750    # Plot if requested 
    758751    view = opts['view'] 
    759752    import matplotlib.pyplot as plt 
    760     if limits is None and not use_data: 
     753    if limits is None: 
    761754        vmin, vmax = np.Inf, -np.Inf 
    762755        if have_base: 
    … …  
    770763    if have_base: 
    771764        if have_comp: plt.subplot(131) 
    772         plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     765        plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
    773766        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    774767        #cbar_title = "log I" 
    … …  
    776769        if have_base: plt.subplot(132) 
    777770        if not opts['is2d'] and have_base: 
    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) 
     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) 
    780773        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    781774        #cbar_title = "log I" 
    … …  
    791784            err[err>cutoff] = cutoff 
    792785        #err,errstr = base/comp,"ratio" 
    793         plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     786        plot_theory(data, None, resid=err, view=errview, use_data=False) 
    794787        if view == 'linear': 
    795788            plt.xscale('linear') 
    … …  
    836829    'linear', 'log', 'q4', 
    837830    'hist', 'nohist', 
    838     'edit', 'html', 'help', 
     831    'edit', 'html', 
    839832    'demo', 'default', 
    840833    ]) 
    841834VALUE_OPTIONS = [ 
    842835    # Note: random is both a name option and a value option 
    843     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 
     836    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 
    844837    ] 
    845838 
    … …  
    947940        'cutoff'    : 0.0, 
    948941        'seed'      : -1,  # default to preset 
    949         'mono'      : True, 
     942        'mono'      : False, 
    950943        # Default to magnetic a magnetic moment is set on the command line 
    951944        'magnetic'  : False, 
    … …  
    958951        'html'      : False, 
    959952        'title'     : None, 
    960         'datafile'  : None, 
    961953    } 
    962954    engines = [] 
    … …  
    979971        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
    980972        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    981         elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    982         elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
     973        elif arg.startswith('-title'):     opts['title'] = arg[7:] 
    983974        elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
    984975        elif arg == '-preset':  opts['seed'] = -1 
    … …  
    1005996        elif arg == '-default':    opts['use_demo'] = False 
    1006997        elif arg == '-html':    opts['html'] = True 
    1007         elif arg == '-help':    opts['html'] = True 
    1008998    # pylint: enable=bad-whitespace 
    1009999 
    … …  
    11221112 
    11231113    # Create the computational engines 
    1124     if opts['datafile'] is not None: 
    1125         data = load_data(os.path.expanduser(opts['datafile'])) 
    1126     else: 
    1127         data, _ = make_data(opts) 
     1114    data, _ = make_data(opts) 
    11281115    if n1: 
    11291116        base = make_engine(model_info, data, engines[0], opts['cutoff']) 
    … …  
    11551142    show html docs for the model 
    11561143    """ 
    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) 
     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 
    11661150 
    11671151def explore(opts): 

  • sasmodels/core.py

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

  • sasmodels/data.py

    r630156b r40a87fa  
    5151    from sas.sascalc.dataloader.loader import Loader  # type: ignore 
    5252    loader = Loader() 
    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: 
     53    data = loader.load(filename) 
     54    if data is None: 
    6155        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 
    7356    return data 
    7457 
    … …  
    365348    # data, but they already handle the masking and graph markup already, so 
    366349    # do not repeat. 
    367     if hasattr(data, 'isSesans') and data.isSesans: 
     350    if hasattr(data, 'lam'): 
    368351        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    369352    elif hasattr(data, 'qx_data'): 
    … …  
    393376    *Iq_calc* is the raw theory values without resolution smearing 
    394377    """ 
    395     if hasattr(data, 'isSesans') and data.isSesans: 
     378    if hasattr(data, 'lam'): 
    396379        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    397380    elif hasattr(data, 'qx_data'): 
    … …  
    463446            if view is 'log': 
    464447                mtheory[mtheory <= 0] = masked 
    465             plt.plot(data.x, scale*mtheory, '-') 
     448            plt.plot(data.x, scale*mtheory, '-', hold=True) 
    466449            all_positive = all_positive and (mtheory > 0).all() 
    467450            some_present = some_present or (mtheory.count() > 0) 
    … …  
    470453            plt.ylim(*limits) 
    471454 
    472         plt.xscale('linear' if not some_present or non_positive_x 
    473                    else view if view is not None 
    474                    else 'log') 
     455        plt.xscale('linear' if not some_present or non_positive_x  else view) 
    475456        plt.yscale('linear' 
    476457                   if view == 'q4' or not some_present or not all_positive 
    477                    else view if view is not None 
    478                    else 'log') 
     458                   else view) 
    479459        plt.xlabel("$q$/A$^{-1}$") 
    480460        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) 
    485461 
    486462    if use_calc: 
    … …  
    502478        if num_plots > 1: 
    503479            plt.subplot(1, num_plots, use_calc + 2) 
    504         plt.plot(data.x, mresid, '.') 
     480        plt.plot(data.x, mresid, '-') 
    505481        plt.xlabel("$q$/A$^{-1}$") 
    506482        plt.ylabel('residuals') 
    507         plt.xscale('linear') 
    508         plt.title('(model - Iq)/dIq') 
     483        plt.xscale('linear' if not some_present or non_positive_x else view) 
    509484 
    510485 
    … …  
    533508        if theory is not None: 
    534509            if is_tof: 
    535                 plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-') 
     510                plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-', hold=True) 
    536511            else: 
    537                 plt.plot(data.x, theory, '-') 
     512                plt.plot(data.x, theory, '-', hold=True) 
    538513        if limits is not None: 
    539514            plt.ylim(*limits) 

  • sasmodels/direct_model.py

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

  • sasmodels/generate.py

    rbb4b509 r1e7b0db0  
    492492 
    493493def test_tag_float(): 
    494     """check that floating point constants are properly identified and tagged with 'f'""" 
    495  
    496     cases = """ 
     494 
     495    cases=""" 
    497496ZP  : 0. 
    498497ZPF : 0.0,0.01,0.1 
    … …  
    520519""" 
    521520 
    522     output = """ 
     521    output=""" 
    523522ZP  : 0.f 
    524523ZPF : 0.0f,0.01f,0.1f 
    … …  
    612611    # type: (str, List[Parameter]) -> List[str] 
    613612    """ 
    614     Return a list of *prefix+parameter* from parameter items. 
    615  
    616     *prefix* should be "v." if v is a struct. 
     613    Return a list of *prefix.parameter* from parameter items. 
    617614    """ 
    618615    return [p.as_call_reference(prefix) for p in pars] 
    … …  
    736733        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
    737734 
    738     magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
     735    magpars = [k-2 for k,p in enumerate(partable.call_parameters) 
    739736               if p.type == 'sld'] 
    740737 
    … …  
    745742    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    746743    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    747     for k, v in enumerate(magpars[:3]): 
     744    for k,v in enumerate(magpars[:3]): 
    748745        source.append("#define MAGNETIC_PAR%d %d"%(k+1, v)) 
    749746 
    … …  
    782779        "#undef CALL_IQ", 
    783780        "#undef KERNEL_NAME", 
    784     ] 
     781         ] 
    785782 
    786783    imagnetic = [ 
    … …  
    875872 
    876873# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog 
    877 RST_PROLOG = r"""\ 
     874RST_PROLOG = """\ 
    878875.. |Ang| unicode:: U+212B 
    879876.. |Ang^-1| replace:: |Ang|\ :sup:`-1` 
    … …  
    931928    from . import rst2html 
    932929    url = "file://"+dirname(info.filename)+"/" 
    933     rst2html.view_html(make_html(info), url=url) 
     930    rst2html.wxview(make_html(info), url=url) 
    934931 
    935932def demo_time(): 

  • sasmodels/kernel_header.c

    rbb4b509 r1e7b0db0  
    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(); on the verified 
     16   // Intel CPU on Mac gives strange values for erf(); also on the tested 
    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: check isnan is correct 
     59         // TODO: test isnan 
    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  
    161150#if 1 
    162151//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
    … …  
    177166#endif 
    178167 
    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 
    199168#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
    200169    q = sqrt(qx*qx + qy*qy); \ 
    … …  
    211180    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
    212181    } while (0) 
    213 #endif 

  • sasmodels/kerneldll.py

    rb2f1d2f r886fa25  
    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 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" 
     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" 
    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

    rdbfd471 rb397165  
    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 = np.asarray(form(), 'd') 
     237            Iq = form() 
    238238            if np.isnan(Iq).any(): continue 
    239239 

  • sasmodels/model_test.py

    rbb4b509 r598857b  
    8585        skip = [] 
    8686    for model_name in models: 
    87         if model_name in skip: 
    88             continue 
     87        if model_name in skip: continue 
    8988        model_info = load_model_info(model_name) 
    9089 
    … …  
    240239            multiple = [test for test in self.info.tests 
    241240                        if isinstance(test[2], list) 
    242                         and not all(result is None for result in test[2])] 
     241                            and not all(result is None for result in test[2])] 
    243242            tests_has_1D_multiple = any(isinstance(test[1][0], float) 
    244243                                        for test in multiple) 
    … …  
    263262            user_pars, x, y = test 
    264263            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)) 
    268264 
    269265            if not isinstance(y, list): 
    … …  
    309305 
    310306    return ModelTestCase 
    311  
    312 def 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  
    328307 
    329308def is_near(target, actual, digits=5): 

  • sasmodels/modelinfo.py

    rbb4b509 rf88e248  
    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         """ 
    348344        # Note: if the parameter is a struct type, then we will need to use 
    349345        # &prefix+id.  For scalars and vectors we can just use prefix+id. 
    … …  
    424420        self.npars = sum(p.length for p in self.kernel_parameters) 
    425421        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
    426                              if p.type == 'sld') 
     422                             if p.type=='sld') 
    427423        self.nvalues = 2 + self.npars 
    428424        if self.nmagnetic: 
    … …  
    461457        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    462458                          for p in self.kernel_parameters) 
    463         self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
     459        self.magnetism_index = [k for k,p in enumerate(self.call_parameters) 
    464460                                if p.id.startswith('M0:')] 
    465461 
    … …  
    548544                              'magnetic', 'magnetic amplitude for '+p.description), 
    549545                    Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], 
    550                               'magnetic', 'magnetic latitude for '+p.description), 
     546                               'magnetic', 'magnetic latitude for '+p.description), 
    551547                    Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], 
    552                               'magnetic', 'magnetic longitude for '+p.description), 
     548                               'magnetic', 'magnetic longitude for '+p.description), 
    553549                ]) 
    554550        #print("call parameters", full_list) 
    … …  
    687683 
    688684    if (model_info.Iq is None 
    689             and model_info.Iqxy is None 
    690             and model_info.Imagnetic is None 
    691             and model_info.form_volume is None): 
     685        and model_info.Iqxy is None 
     686        and model_info.Imagnetic is None 
     687        and model_info.form_volume is None): 
    692688        return 
    693689 
    … …  
    847843    #: it to False because they require double precision calculations. 
    848844    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 
    853845    #: True if the model is a structure factor used to model the interaction 
    854846    #: between form factor models.  This will default to False if it is not 

  • sasmodels/models/barbell.py

    r9802ab3 rfcb33e4  
    6868The 2D scattering intensity is calculated similar to the 2D cylinder model. 
    6969 
    70 .. figure:: img/cylinder_angle_definition.png 
     70.. figure:: img/cylinder_angle_definition.jpg 
    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, sin, cos, pi 
     89from numpy import inf 
    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, [-360, 360], "orientation", "Barbell axis to beam angle"], 
    111               ["phi",         "degrees",     60, [-360, 360], "orientation", "Rotation about beam"], 
     110              ["theta",       "degrees",     60, [-inf, inf], "orientation", "In plane angle"], 
     111              ["phi",         "degrees",     60, [-inf, inf], "orientation", "Out of plane angle"], 
    112112             ] 
    113113# pylint: enable=bad-whitespace, line-too-long 
    … …  
    125125            phi_pd=15, phi_pd_n=0, 
    126126           ) 
    127 q = 0.1 
    128 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    129 qx = q*cos(pi/6.0) 
    130 qy = q*sin(pi/6.0) 
    131 tests = [[{}, 0.075, 25.5691260532], 
    132         [{'theta':80., 'phi':10.}, (qx, qy), 3.04233067789], 
    133         ] 
    134 del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/bcc_paracrystal.c

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

  • sasmodels/models/bcc_paracrystal.py

    r69e1afc r925ad6e  
    7979be accurate. 
    8080 
    81 .. figure:: img/parallelepiped_angle_definition.png 
     81.. figure:: img/bcc_angle_definition.png 
    8282 
    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). 
     83    Orientation of the crystal with respect to the scattering plane. 
    8584 
    8685References 
    … …  
    10099""" 
    101100 
    102 from numpy import inf, pi 
     101from numpy import inf 
    103102 
    104103name = "bcc_paracrystal" 
    … …  
    123122              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "sld",         "Particle scattering length density"], 
    124123              ["sld_solvent", "1e-6/Ang^2",  1,    [-inf, inf], "sld",         "Solvent scattering length density"], 
    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"] 
     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"] 
    128127             ] 
    129128# pylint: enable=bad-whitespace, line-too-long 
    … …  
    142141    psi_pd=15, psi_pd_n=0, 
    143142    ) 
    144 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    145 # add 2d test later 
    146 q =4.*pi/220. 
    147 tests = [ 
    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

    r9802ab3 rfcb33e4  
    7171The 2D scattering intensity is calculated similar to the 2D cylinder model. 
    7272 
    73 .. figure:: img/cylinder_angle_definition.png 
     73.. figure:: img/cylinder_angle_definition.jpg 
    7474 
    7575    Definition of the angles for oriented 2D cylinders. 
    … …  
    9191 
    9292""" 
    93 from numpy import inf, sin, cos, pi 
     93from numpy import inf 
    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, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    132               ["phi",        "degrees", 60, [-360, 360], "orientation", "rotation about beam"], 
     131              ["theta",      "degrees", 60, [-inf, inf], "orientation", "inclination angle"], 
     132              ["phi",        "degrees", 60, [-inf, inf], "orientation", "deflection angle"], 
    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) 
    147 q = 0.1 
    148 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    149 qx = q*cos(pi/6.0) 
    150 qy = q*sin(pi/6.0) 
    151 tests = [[{}, 0.075, 26.0698570695], 
    152         [{'theta':80., 'phi':10.}, (qx, qy), 0.561811990502], 
    153         ] 
    154 del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/core_shell_bicelle.c

    rb260926 r592343f  
    3030 
    3131static double 
    32 bicelle_kernel(double q, 
     32bicelle_kernel(double qq, 
    3333              double rad, 
    3434              double radthick, 
    3535              double facthick, 
    36               double halflength, 
     36              double length, 
    3737              double rhoc, 
    3838              double rhoh, 
    … …  
    4242              double cos_alpha) 
    4343{ 
     44    double si1,si2,be1,be2; 
     45 
    4446    const double dr1 = rhoc-rhoh; 
    4547    const double dr2 = rhor-rhosolv; 
    4648    const double dr3 = rhoh-rhor; 
    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); 
     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; 
    5056 
    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); 
     57    be1 = sas_2J1x_x(besarg1); 
     58    be2 = sas_2J1x_x(besarg2); 
     59    si1 = sas_sinx_x(sinarg1); 
     60    si2 = sas_sinx_x(sinarg2); 
    5561 
    5662    const double t = vol1*dr1*si1*be1 + 
    … …  
    5864                     vol3*dr3*si2*be1; 
    5965 
    60     const double retval = t*t; 
     66    const double retval = t*t*sin_alpha; 
    6167 
    6268    return retval; 
    … …  
    6571 
    6672static double 
    67 bicelle_integration(double q, 
     73bicelle_integration(double qq, 
    6874                   double rad, 
    6975                   double radthick, 
    … …  
    7783    // set up the integration end points 
    7884    const double uplim = M_PI_4; 
    79     const double halflength = 0.5*length; 
     85    const double halfheight = 0.5*length; 
    8086 
    8187    double summ = 0.0; 
    … …  
    8490        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8591        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
    87                              halflength, rhoc, rhoh, rhor, rhosolv, 
     92        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
     93                             halfheight, rhoc, rhoh, rhor, rhosolv, 
    8894                             sin_alpha, cos_alpha); 
    89         summ += yyy*sin_alpha; 
     95        summ += yyy; 
    9096    } 
    9197 
    … …  
    113119    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
    114120                           0.5*length, core_sld, face_sld, rim_sld, 
    115                            solvent_sld, sin_alpha, cos_alpha); 
    116     return 1.0e-4*answer; 
     121                           solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 
     122 
     123    answer *= 1.0e-4; 
     124 
     125    return answer; 
    117126} 
    118127 

  • sasmodels/models/core_shell_bicelle.py

    r9802ab3 r8afefae  
    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 For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data,  
    66 see the :ref:`cylinder` model for further information. 
     65The *theta* and *phi* parameters are not used for the 1D output. 
    6766Our implementation of the scattering kernel and the 1D scattering intensity 
    6867use the c-library from NIST. 
    6968 
    70 .. figure:: img/cylinder_angle_definition.png 
     69.. figure:: img/cylinder_angle_definition.jpg 
    7170 
    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$. 
     71    Definition of the angles for the oriented core shell bicelle tmodel. 
    7572 
    7673 
    … …  
    9188""" 
    9289 
    93 from numpy import inf, sin, cos, pi 
     90from numpy import inf, sin, cos 
    9491 
    9592name = "core_shell_bicelle" 
    … …  
    138135    ["sld_rim",        "1e-6/Ang^2", 4, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    139136    ["sld_solvent",    "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Solvent scattering length density"], 
    140     ["theta",          "degrees",   90, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    141     ["phi",            "degrees",    0, [-360, 360], "orientation", "rotation about beam"] 
     137    ["theta",          "degrees",   90, [-inf, inf], "orientation", "In plane angle"], 
     138    ["phi",            "degrees",    0, [-inf, inf], "orientation", "Out of plane angle"], 
    142139    ] 
    143140 
    … …  
    158155            theta=90, 
    159156            phi=0) 
    160 q = 0.1 
    161 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    162 qx = q*cos(pi/6.0) 
    163 qy = q*sin(pi/6.0) 
    164 tests = [[{}, 0.05, 7.4883545957], 
    165         [{'theta':80., 'phi':10.}, (qx, qy), 2.81048892474 ] 
    166         ] 
    167 del qx, qy  # not necessary to delete, but cleaner 
    168157 
     158#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rdedcf34 r592343f  
     1double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length); 
     2double 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 
     14double 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 
    128// NOTE that "length" here is the full height of the core! 
    2 static double 
    3 form_volume(double r_minor, 
    4         double x_core, 
    5         double thick_rim, 
    6         double thick_face, 
    7         double length) 
     29double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length) 
    830{ 
    9     return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 
     31    return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face); 
    1032} 
    1133 
    12 static double 
    13 Iq(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) 
     34double  
     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) 
    2345{ 
    2446    double si1,si2,be1,be2; 
    2547     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    26      // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
     48     // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle 
    2749     //    const double uplim = M_PI_4; 
    2850    const double halfheight = 0.5*length; 
    … …  
    3355    //const double vbj=M_PI; 
    3456 
    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); 
     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); 
    4466 
    4567    //initialize integral 
    … …  
    5274        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    5375        double inner_sum=0; 
    54         double sinarg1 = q*halfheight*cos_alpha; 
    55         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
     76        double sinarg1 = qq*halfheight*cos_alpha; 
     77        double sinarg2 = qq*(halfheight+facthick)*cos_alpha; 
    5678        si1 = sas_sinx_x(sinarg1); 
    5779        si2 = sas_sinx_x(sinarg2); 
    … …  
    6082            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    6183            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    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; 
     84            const double rr = sqrt(rA - rB*cos(beta)); 
     85            double besarg1 = qq*rr*sin_alpha; 
     86            double besarg2 = qq*(rr+radthick)*sin_alpha; 
    6587            be1 = sas_2J1x_x(besarg1); 
    6688            be2 = sas_2J1x_x(besarg2); 
    … …  
    7698} 
    7799 
    78 static double 
     100double  
    79101Iqxy(double qx, double qy, 
    80           double r_minor, 
     102          double rad, 
    81103          double x_core, 
    82           double thick_rim, 
    83           double thick_face, 
     104          double radthick, 
     105          double facthick, 
    84106          double length, 
    85107          double rhoc, 
    … …  
    92114{ 
    93115       // THIS NEEDS TESTING 
    94     double q, xhat, yhat, zhat; 
    95     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     116    double qq, cos_val, cos_mu, cos_nu; 
     117    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu); 
    96118    const double dr1 = rhoc-rhoh; 
    97119    const double dr2 = rhor-rhosolv; 
    98120    const double dr3 = rhoh-rhor; 
    99     const double r_major = r_minor*x_core; 
     121    const double radius_major = rad*x_core; 
    100122    const double halfheight = 0.5*length; 
    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); 
     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); 
    104126 
    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 ); 
     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 ); 
    113135    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); 
    114137    return 1.0e-4 * Aq; 
    115138} 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r9802ab3 r8afefae  
    7676bicelles is then given by integrating over all possible $\alpha$ and $\psi$. 
    7777 
    78 For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     78For oriented bicellles the *theta*, *phi* and *psi* orientation parameters only appear when fitting 2D data,  
    7979see the :ref:`elliptical-cylinder` model for further information. 
    8080 
    8181 
    82 .. figure:: img/elliptical_cylinder_angle_definition.png 
     82.. figure:: img/elliptical_cylinder_angle_definition.jpg 
    8383 
    84     Definition of the angles for the oriented core_shell_bicelle_elliptical particles.    
    85  
     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. 
    8686 
    8787 
    … …  
    9999""" 
    100100 
    101 from numpy import inf, sin, cos, pi 
     101from numpy import inf, sin, cos 
    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.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"] 
     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"], 
    131131    ] 
    132132 
    … …  
    150150            psi=0) 
    151151 
    152 q = 0.1 
    153 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
    154 qx = q*cos(pi/6.0) 
    155 qy = q*sin(pi/6.0) 
     152#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 
    156153 
    157154tests = [ 
    … …  
    162159    'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 
    163160    0.015, 286.540286], 
    164 #    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
    165         ] 
    166  
    167 del qx, qy  # not necessary to delete, but cleaner 
     161] 

  • sasmodels/models/core_shell_cylinder.py

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

  • sasmodels/models/core_shell_ellipsoid.py

    r9802ab3 r8e68ea0  
    7777   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    7878 
    79 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
    80 see the :ref:`elliptical-cylinder` model for further information. 
    8179 
    8280References 
    … …  
    134132    ["sld_shell",     "1e-6/Ang^2", 1,   [-inf, inf], "sld",         "Shell scattering length density"], 
    135133    ["sld_solvent",   "1e-6/Ang^2", 6.3, [-inf, inf], "sld",         "Solvent scattering length density"], 
    136     ["theta",         "degrees",    0,   [-360, 360], "orientation", "elipsoid axis to beam angle"], 
    137     ["phi",           "degrees",    0,   [-360, 360], "orientation", "rotation about beam"], 
     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"], 
    138136    ] 
    139137# pylint: enable=bad-whitespace, line-too-long 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c r1e7b0db0  
    134134    double psi) 
    135135{ 
    136     double q, zhat, yhat, xhat; 
    137     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     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); 
    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*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); 
     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); 
    168168     
    169169 

  • sasmodels/models/core_shell_parallelepiped.py

    r9b79f29 rcb0dc22  
    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 1D calculation does **NOT** actually calculate a c face rim despite the presence of 
    9 the parameter. Some other aspects of the 1D calculation may be wrong. 
     8the model does **NOT** actually calculate a c face rim despite the presence of 
     9the parameter. 
    1010 
    1111.. note:: 
    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. 
     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. 
    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] 
    5750 
    5851.. note:: 
    5952 
    60     Why does t_B not appear in the above equation? 
    6153    For the calculation of the form factor to be valid, the sides of the solid 
    62     MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
     54    MUST be chosen such that** $A < B < C$. 
    6355    If this inequality is not satisfied, the model will not report an error, 
    6456    but the calculation will not be correct and thus the result wrong. 
    6557 
    6658FITTING NOTES 
    67 If the scale is set equal to the particle volume fraction, $\phi$, the returned 
     59If the scale is set equal to the particle volume fraction, |phi|, the returned 
    6860value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    6961However, **no interparticle interference effects are included in this 
    … …  
    8173NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    8274based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    83 and length $(C+2t_C)$ values, after appropriately 
    84 sorting the three dimensions to give an oblate or prolate particle, to give an  
    85 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     75and length $(C+2t_C)$ values, and used as the effective radius 
     76for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    8677 
    8778To provide easy access to the orientation of the parallelepiped, we define the 
    … …  
    9283*x*-axis of the detector. 
    9384 
    94 .. figure:: img/parallelepiped_angle_definition.png 
     85.. figure:: img/parallelepiped_angle_definition.jpg 
    9586 
    9687    Definition of the angles for oriented core-shell parallelepipeds. 
    9788 
    98 .. figure:: img/parallelepiped_angle_projection.png 
     89.. figure:: img/parallelepiped_angle_projection.jpg 
    9990 
    10091    Examples of the angles for oriented core-shell parallelepipeds against the 
    … …  
    121112 
    122113import numpy as np 
    123 from numpy import pi, inf, sqrt, cos, sin 
     114from numpy import pi, inf, sqrt 
    124115 
    125116name = "core_shell_parallelepiped" 
    … …  
    153144              ["thick_rim_c", "Ang", 10, [0, inf], "volume", 
    154145               "Thickness of C rim"], 
    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"], 
     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"], 
    161152             ] 
    162153 
    … …  
    195186            psi_pd=10, psi_pd_n=1) 
    196187 
    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 
    198 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     188qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    199189tests = [[{}, 0.2, 0.533149288477], 
    200190         [{}, [0.2], [0.533149288477]], 
    201          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
    202          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
     191         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.032102135569], 
     192         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.032102135569]], 
    203193        ] 
    204194del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/cylinder.py

    r9802ab3 rb7e8b94  
    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.png 
     63.. figure:: img/cylinder_angle_definition.jpg 
    6464 
    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. 
     65    Definition of the angles for oriented cylinders. 
    6966 
    70 .. figure:: img/cylinder_angle_projection.png 
    71  
    72     Examples for oriented cylinders. 
    73  
    74 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.  
    75 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel  
    77 to 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  
    79 situations, with very broad distributions, should still be approached with care.)  
     67The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
    8068 
    8169Validation 
    … …  
    135123              ["length", "Ang", 400, [0, inf], "volume", 
    136124               "Cylinder length"], 
    137               ["theta", "degrees", 60, [-360, 360], "orientation", 
    138                "cylinder axis to beam angle"], 
    139               ["phi", "degrees",   60, [-360, 360], "orientation", 
    140                "rotation about beam"], 
     125              ["theta", "degrees", 60, [-inf, inf], "orientation", 
     126               "latitude"], 
     127              ["phi", "degrees", 60, [-inf, inf], "orientation", 
     128               "longitude"], 
    141129             ] 
    142130 
    … …  
    164152tests = [[{}, 0.2, 0.042761386790780453], 
    165153        [{}, [0.2], [0.042761386790780453]], 
    166 #  new coords 
     154#  new coords     
    167155        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
    168156        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 

  • sasmodels/models/ellipsoid.c

    r3b571ae r130d4c7  
    33double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
     5 
     6static 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} 
    527 
    628double form_volume(double radius_polar, double radius_equatorial) 
    … …  
    1537    double radius_equatorial) 
    1638{ 
    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  
    2639    // translate a point in [-1,1] to a point in [0, 1] 
    27     // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2840    const double zm = 0.5; 
    2941    const double zb = 0.5; 
    3042    double total = 0.0; 
    3143    for (int i=0;i<76;i++) { 
    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; 
     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); 
    3647    } 
    3748    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5162    double q, sin_alpha, cos_alpha; 
    5263    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, 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); 
     64    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    5665    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5766 
    58     return 1.0e-4 * square(f * s); 
     67    return 1.0e-4 * form * s * s; 
    5968} 
    6069 

  • sasmodels/models/ellipsoid.py

    r9b79f29 r925ad6e  
    1818.. math:: 
    1919 
    20     F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3} 
     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} 
    2123 
    22 for 
     24and 
    2325 
    2426.. math:: 
    2527 
    26     r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2} 
     28    r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha 
     29        + R_p^2 \cos^2 \alpha \right]^{1/2} 
    2730 
    2831 
    2932$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$, 
    30 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar 
    31 radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 
    32 radius perpendicular to the rotational axis of the ellipsoid and 
    33 $\Delta \rho$ (contrast) is the scattering length density difference between 
    34 the scatterer and the solvent. 
     33$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the 
     34rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular 
     35to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the 
     36scattering length density difference between the scatterer and the solvent. 
    3537 
    36 For randomly oriented particles use the orientational average, 
     38For randomly oriented particles: 
    3739 
    3840.. math:: 
    3941 
    40    \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha} 
     42   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    4143 
    42  
    43 computed 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  
    49 with 
    50  
    51 .. math:: 
    52  
    53     r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2} 
    5444 
    5545To provide easy access to the orientation of the ellipsoid, we define 
    … …  
    5848:ref:`cylinder orientation figure <cylinder-angle-definition>`. 
    5949For the ellipsoid, $\theta$ is the angle between the rotational axis 
    60 and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 
    61 in the $xy$ plane. 
     50and the $z$ -axis. 
    6251 
    6352NB: The 2nd virial coefficient of the solid ellipsoid is calculated based 
    … …  
    10190than 500. 
    10291 
    103 Model was also tested against the triaxial ellipsoid model with equal major 
    104 and minor equatorial radii.  It is also consistent with the cyclinder model 
    105 with polar radius equal to length and equatorial radius equal to radius. 
    106  
    10792References 
    10893---------- 
    … …  
    11196*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
    11297Plenum Press, New York, 1987. 
    113  
    114 Authorship 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 
    12098""" 
    12199 
    122 from numpy import inf, sin, cos, pi 
     100from numpy import inf 
    123101 
    124102name = "ellipsoid" 
    … …  
    151129              ["radius_equatorial", "Ang", 400, [0, inf], "volume", 
    152130               "Equatorial radius"], 
    153               ["theta", "degrees", 60, [-360, 360], "orientation", 
    154                "ellipsoid axis to beam angle"], 
    155               ["phi", "degrees", 60, [-360, 360], "orientation", 
    156                "rotation about beam"], 
     131              ["theta", "degrees", 60, [-inf, inf], "orientation", 
     132               "In plane angle"], 
     133              ["phi", "degrees", 60, [-inf, inf], "orientation", 
     134               "Out of plane angle"], 
    157135             ] 
    158136 
    … …  
    161139def ER(radius_polar, radius_equatorial): 
    162140    import numpy as np 
    163 # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     141 
    164142    ee = np.empty_like(radius_polar) 
    165143    idx = radius_polar > radius_equatorial 
    … …  
    190168            theta_pd=15, theta_pd_n=45, 
    191169            phi_pd=15, phi_pd_n=1) 
    192 q = 0.1 
    193 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    194 qx = q*cos(pi/6.0) 
    195 qy = q*sin(pi/6.0) 
    196 tests = [[{}, 0.05, 54.8525847025], 
    197         [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 ], 
    198         ] 
    199 del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/elliptical_cylinder.c

    r61104c8 r592343f  
    6767     double theta, double phi, double psi) 
    6868{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     69    double q, cos_val, cos_mu, cos_nu; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 
    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*xhat) + square(yhat)); 
     74    const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 
    7575    const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     76    const double si = sas_sinx_x(q*0.5*length*cos_val); 
    7777    const double Aq = be * si; 
    7878    const double delrho = sld - solvent_sld; 

  • sasmodels/models/elliptical_cylinder.py

    r9802ab3 rfcb33e4  
    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.  
     59$\Psi$ is the rotational angle around its own long_c axis against the $q$ plane. 
     60For example, $\Psi = 0$ when the $r_\text{minor}$ axis is parallel to the 
     61$x$ axis of the detector. 
    6062 
    6163All angle parameters are valid and given only for 2D calculation; ie, an 
    6264oriented system. 
    6365 
    64 .. figure:: img/elliptical_cylinder_angle_definition.png 
     66.. figure:: img/elliptical_cylinder_angle_definition.jpg 
    6567 
    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. 
     68    Definition of angles for 2D 
    6869 
    69 .. figure:: img/elliptical_cylinder_angle_projection.png 
     70.. figure:: img/cylinder_angle_projection.jpg 
    7071 
    7172    Examples of the angles for oriented elliptical cylinders against the 
    72     detector plane, with $\Psi$ = 0. 
    73  
    74 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.  
    75 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the  
    77 cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)  
    78 The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to  
    79 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation  
    80 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)  
     73    detector plane. 
    8174 
    8275NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    … …  
    115108""" 
    116109 
    117 from numpy import pi, inf, sqrt, sin, cos 
     110from numpy import pi, inf, sqrt 
    118111 
    119112name = "elliptical_cylinder" 
    … …  
    132125              ["sld",         "1e-6/Ang^2", 4.0,   [-inf, inf], "sld",         "Cylinder scattering length density"], 
    133126              ["sld_solvent", "1e-6/Ang^2", 1.0,   [-inf, inf], "sld",         "Solvent scattering length density"], 
    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"]] 
     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"]] 
    137130 
    138131# pylint: enable=bad-whitespace, line-too-long 
    … …  
    156149                           + (length + radius) * (length + pi * radius)) 
    157150    return 0.5 * (ddd) ** (1. / 3.) 
    158 q = 0.1 
    159 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
    160 qx = q*cos(pi/6.0) 
    161 qy = q*sin(pi/6.0) 
    162151 
    163152tests = [ 
    … …  
    169158      'sld_solvent':1.0, 'background':0.0}, 
    170159     0.001, 675.504402], 
    171 #    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
    172160] 

  • sasmodels/models/fcc_paracrystal.c

    r50beefe r4962519  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    9494 
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
     95    const double a1 = cos_a2 + cos_a3; 
     96    const double a2 = cos_a3 + cos_a1; 
     97    const double a3 = cos_a2 + cos_a1; 
    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

    r69e1afc r925ad6e  
    76762D model computation. 
    7777 
    78 .. figure:: img/parallelepiped_angle_definition.png 
     78.. figure:: img/bcc_angle_definition.png 
    7979 
    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). 
     80    Orientation of the crystal with respect to the scattering plane. 
    8281 
    8382References 
    … …  
    9190""" 
    9291 
    93 from numpy import inf, pi 
     92from numpy import inf 
    9493 
    9594name = "fcc_paracrystal" 
    … …  
    111110              ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], 
    112111              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"], 
    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"] 
     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"] 
    116115             ] 
    117116# pylint: enable=bad-whitespace, line-too-long 
    … …  
    129128            psi_pd=15, psi_pd_n=0, 
    130129           ) 
    131 # april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    132 q =4.*pi/220. 
    133 tests = [ 
    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

    r5702f52 ra807206  
    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, 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"], 
     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"], 
    7979    ["dielectconst",  "None",    71.08,   [-inf, inf], "", "dielectric constant (relative permittivity) of solvent, default water, for Debye length"] 
    8080    ] 

  • sasmodels/models/hollow_cylinder.py

    r9b79f29 raea2e2a  
    6060""" 
    6161 
    62 from numpy import pi, inf, sin, cos 
     62from numpy import pi, inf 
    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", "Cylinder axis to beam angle"], 
    85     ["phi",         "degrees",  0,   [-360, 360], "orientation", "Rotation about beam"], 
     84    ["theta",       "degrees", 90,   [-360, 360], "orientation", "Theta angle"], 
     85    ["phi",         "degrees",  0,   [-360, 360], "orientation", "Phi angle"], 
    8686    ] 
    8787# pylint: enable=bad-whitespace, line-too-long 
    … …  
    129129            theta_pd=10, theta_pd_n=5, 
    130130           ) 
    131 q = 0.1 
    132 # april 6 2017, rkh added a 2d unit test, assume correct! 
    133 qx = q*cos(pi/6.0) 
    134 qy = q*sin(pi/6.0) 
     131 
    135132# Parameters for unit tests 
    136133tests = [ 
    137134    [{}, 0.00005, 1764.926], 
    138135    [{}, 'VR', 1.8], 
    139     [{}, 0.001, 1756.76], 
    140     [{}, (qx, qy), 2.36885476192  ], 
    141         ] 
    142 del qx, qy  # not necessary to delete, but cleaner 
     136    [{}, 0.001, 1756.76] 
     137    ] 

  • sasmodels/models/lib/sas_3j1x_x.c

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

  • sasmodels/models/lib/sas_J0.c

    reb2946f rc8902ac  
    236236        xx = x; 
    237237 
    238     // 2017-05-18 PAK - support negative x 
    239     if( xx <= 2.0 ) { 
     238    if( x <= 2.0 ) { 
    240239        z = xx * xx; 
    241         if( xx < 1.0e-3 ) 
     240        if( x < 1.0e-3 ) 
    242241            return( 1.0 - 0.25*z ); 
    243242 
    … …  
    246245    } 
    247246 
    248     q = 1.0/xx; 
     247    q = 1.0/x; 
    249248    w = sqrt(q); 
    250249 

  • sasmodels/models/lib/sas_J1.c

    reb2946f r473a9f1  
    109109{ 
    110110 
    111     double w, z, p, q, xn, abs_x, sign_x; 
     111    double w, z, p, q, xn; 
    112112 
    113113    const double Z1 = 1.46819706421238932572E1; 
    … …  
    116116    const double SQ2OPI = 0.79788456080286535588; 
    117117 
    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; 
     118    w = x; 
     119    if( x < 0 ) 
     120        w = -x; 
     121 
     122    if( w <= 5.0 ) { 
     123        z = x * x; 
     124        w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
     125        w = w * x * (z - Z1) * (z - Z2); 
     126        return( w ); 
    125127    } 
    126128 
    127     if( abs_x <= 5.0 ) { 
    128         z = abs_x * abs_x; 
    129         w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
    130         w = w * abs_x * (z - Z1) * (z - Z2); 
    131         return( sign_x * w ); 
    132     } 
    133  
    134     w = 5.0/abs_x; 
     129    w = 5.0/x; 
    135130    z = w * w; 
     131 
    136132    p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 
    137133    q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 
    138     xn = abs_x - THPIO4; 
     134 
     135    xn = x - THPIO4; 
    139136 
    140137    double sn, cn; 
    … …  
    142139    p = p * cn - w * q * sn; 
    143140 
    144     return( sign_x * p * SQ2OPI / sqrt(abs_x) ); 
     141    return( p * SQ2OPI / sqrt(x) ); 
    145142} 
    146143 
    … …  
    182179    }; 
    183180 
    184 float cephes_j1f(float xx) 
     181float cephes_j1f(float x) 
    185182{ 
    186183 
    187     float x, w, z, p, q, xn; 
     184    float xx, w, z, p, q, xn; 
    188185 
    189186    const float Z1 = 1.46819706421238932572E1; 
    … …  
    191188 
    192189 
    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 ); 
     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 ); 
    202198    } 
    203199 
    … …  
    208204    w = q*q; 
    209205    xn = q * polevl( w, PH1J1, 7) - THPIO4F; 
    210     p = p * cos(xn + x); 
    211  
    212     return( xx < 0. ? -p : p ); 
     206    p = p * cos(xn + xx); 
     207 
     208    return(p); 
    213209} 
    214210#endif 

  • sasmodels/models/lib/sas_erf.c

    reb2946f rb3796fa  
    165165        // the erf function instead and z < 1.0. 
    166166        //return (1.0 - cephes_erf(a)); 
    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; 
     167        z = x * x; 
     168        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5); 
     169 
     170        return y; 
    172171    } 
    173172 
    … …  
    280279        //is explicit here for the case < 1.0 
    281280        //return (1.0 - sas_erf(a)); 
    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; 
     281        z = x * x; 
     282        y = x * polevl( z, TF, 6 ); 
     283 
     284        return y; 
    287285    } 
    288286 

  • sasmodels/models/onion.py

    rbccb40f r925ad6e  
    11r""" 
    22This model provides the form factor, $P(q)$, for a multi-shell sphere where 
    3 the scattering length density (SLD) of each shell is described by an 
     3the scattering length density (SLD) of the 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 in this case) and becomes flat. We set the constant 
     144$\rho_\text{out}$ is ignored 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_in[k]) 
     348            rho.append(sld_out[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

    rd605080 r1e7b0db0  
    6767    double psi) 
    6868{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     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); 
    7171 
    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); 
     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); 
    7575    const double V = form_volume(length_a, length_b, length_c); 
    7676    const double drho = (sld - solvent_sld); 

  • sasmodels/models/parallelepiped.py

    r34a9e4e r8e68ea0  
    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  
    1817.. figure:: img/parallelepiped_geometry.jpg 
    1918 
    … …  
    2221.. note:: 
    2322 
    24 The three dimensions of the parallelepiped (strictly here a cuboid) may be given in  
    25 $any$ size order. To avoid multiple fit solutions, especially 
    26 with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may  
    27 be a number of closely similar "best fits", so some trial and error, or fixing of some  
    28 dimensions at expected values, may help. 
     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. 
    2926 
    3027The 1D scattering intensity $I(q)$ is calculated as: 
    … …  
    7067    \mu &= qB 
    7168 
     69 
    7270The scattering intensity per unit volume is returned in units of |cm^-1|. 
    7371 
    7472NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    75 the averaged effective radius, after appropriately sorting the three 
    76 dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
     73the averaged effective radius $(=\sqrt{A B / \pi})$ and 
    7774length $(= C)$ values, and used as the effective radius for 
    7875$S(q)$ when $P(q) \cdot S(q)$ is applied. 
    … …  
    105102.. _parallelepiped-orientation: 
    106103 
    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 
     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 
    114111    detector plane. 
    115112 
    116 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    117 appear. 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  
    119 about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to  
    120 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation  
    121 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)  
    122  
    123      
    124113For a given orientation of the parallelepiped, the 2D form factor is 
    125114calculated as 
    … …  
    127116.. math:: 
    128117 
    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 
     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 
    132121 
    133122with 
    … …  
    165154angles. 
    166155 
     156This model is based on form factor calculations implemented in a c-library 
     157provided by the NIST Center for Neutron Research (Kline, 2006). 
    167158 
    168159References 
    … …  
    172163 
    173164R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    174  
    175 Authorship 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  
    183165""" 
    184166 
    185167import numpy as np 
    186 from numpy import pi, inf, sqrt, sin, cos 
     168from numpy import pi, inf, sqrt 
    187169 
    188170name = "parallelepiped" 
    … …  
    198180            mu = q*B 
    199181        V: Volume of the rectangular parallelepiped 
    200         alpha: angle between the long axis of the 
     182        alpha: angle between the long axis of the  
    201183            parallelepiped and the q-vector for 1D 
    202184""" 
    … …  
    214196              ["length_c", "Ang", 400, [0, inf], "volume", 
    215197               "Larger side of the parallelepiped"], 
    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"], 
     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"], 
    222204             ] 
    223205 
    … …  
    226208def ER(length_a, length_b, length_c): 
    227209    """ 
    228     Return effective radius (ER) for P(q)*S(q) 
     210        Return effective radius (ER) for P(q)*S(q) 
    229211    """ 
    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]) 
     212 
    236213    # surface average radius (rough approximation) 
    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)) 
     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)) 
    240217    return 0.5 * (ddd) ** (1. / 3.) 
    241218 
    … …  
    253230            phi_pd=10, phi_pd_n=1, 
    254231            psi_pd=10, psi_pd_n=10) 
    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 
    256 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     232 
     233qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    257234tests = [[{}, 0.2, 0.17758004974], 
    258235         [{}, [0.2], [0.17758004974]], 
    259          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475], 
    260          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]], 
     236         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014], 
     237         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]], 
    261238        ] 
    262239del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/sc_paracrystal.c

    r50beefe r4962519  
    111111          double psi) 
    112112{ 
    113     double q, zhat, yhat, xhat; 
    114     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     113    double q, cos_a1, cos_a2, cos_a3; 
     114    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    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*zhat)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
     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); 
    123123 
    124124    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 

  • sasmodels/models/sc_paracrystal.py

    r69e1afc r925ad6e  
    8383model computation. 
    8484 
    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). 
     85.. figure:: img/sc_crystal_angle_definition.jpg 
    8986 
    9087Reference 
    … …  
    130127              ["sld",  "1e-6/Ang^2",         3.0, [0.0, inf],  "sld",         "Sphere scattering length density"], 
    131128              ["sld_solvent", "1e-6/Ang^2",  6.3, [0.0, inf],  "sld",         "Solvent scattering length density"], 
    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"] 
     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"], 
    135132             ] 
    136133# pylint: enable=bad-whitespace, line-too-long 
    … …  
    149146 
    150147tests = [ 
    151     # Accuracy tests based on content in test/utest_extra_models.py, 2d tests added April 10, 2017 
     148    # Accuracy tests based on content in test/utest_extra_models.py 
    152149    [{}, 0.001, 10.3048], 
    153150    [{}, 0.215268, 0.00814889], 
    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 ] 
     151    [{}, (0.414467), 0.001313289] 
    157152    ] 
    158153 

  • sasmodels/models/spherical_sld.py

    r63a7fe8 r925ad6e  
    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

    r1b85b55 rb7e8b94  
    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.png 
     79.. figure:: img/cylinder_angle_definition.jpg 
    8080 
    8181    Examples of the angles against the detector plane. 
    … …  
    103103""" 
    104104 
    105 from numpy import inf, sin, cos, pi 
     105from numpy import inf 
    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,   [-360, 360], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"], 
    134     ["phi",         "degrees",     0,   [-360, 360], "orientation", "Rotation about beam"], 
     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"], 
    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 q = 0.1 
    155 # april 6 2017, rkh added a 2d unit test, assume correct! 
    156 qx = q*cos(pi/6.0) 
    157 qy = q*sin(pi/6.0) 
     154# but should not matter here as so far all the tests are 1D not 2D 
     155tests = [ 
    158156# 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 
    162 tests = [ 
    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, 
     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, 
    171166      'theta': 90.0, 
    172167      'phi': 0.0, 
    … …  
    181176      'sld_core': 4.0, 
    182177      'sld_layer': -0.4, 
    183       'sld_solvent': 5.0, 
     178      'solvent_sd': 5.0, 
    184179      'theta': 90.0, 
    185180      'phi': 0.0, 
    … …  
    191186    [{'thick_core': 10.0, 
    192187      'thick_layer': 15.0, 
    193       'radius': 100.0, 
     188      'radius': 3000.0, 
    194189      'n_stacking': 5, 
    195190      'sigma_d': 0.0, 
    196191      'sld_core': 4.0, 
    197192      'sld_layer': -0.4, 
    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, 
     193      'solvent_sd': 5.0, 
    212194      'theta': 90.0, 
    213195      'phi': 0.0, 
    … …  
    225207      'sld_core': 4.0, 
    226208      'sld_layer': -0.4, 
    227       'sld_solvent': 5.0, 
     209      'solvent_sd': 5.0, 
    228210      'theta': 90.0, 
    229211      'phi': 0.0, 
    … …  
    240222      'sld_core': 4.0, 
    241223      'sld_layer': -0.4, 
    242       'sld_solvent': 5.0, 
     224      'solvent_sd': 5.0, 
    243225      'theta': 90.0, 
    244226      'phi': 0.0, 
    … …  
    246228      'background': 0.001, 
    247229     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
    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, 
     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, 
    272239      'theta': 90.0, 
    273240      'phi': 0.0, 
    … …  
    276243     }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
    277244    ] 
    278 # 11Jan2017   RKH checking unit test again, note they are all 1D, no 2D 
     245# 11Jan2017   RKH checking unit test agai, note they are all 1D, no 2D 
     246 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 r925ad6e  
    2020    double radius_polar) 
    2121{ 
    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; 
     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; 
    2726    double outer = 0.0; 
    2827    for (int i=0;i<76;i++) { 
    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)); 
     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; 
    3234 
    3335        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             // 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; 
     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 ; 
    4241        } 
    43         outer += Gauss76Wt[i] * inner;  // correcting for dx later 
     42        outer += Gauss76Wt[i] * 0.5 * inner; 
    4443    } 
    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); 
     44    // translate dx in [-1,1] to dx in [lower,upper] 
     45    const double fqsq = outer*zm; 
    4746    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    4847    return 1.0e-4 * s * s * fqsq; 
    … …  
    5958    double psi) 
    6059{ 
    61     double q, xhat, yhat, zhat; 
    62     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     60    double q, calpha, cmu, cnu; 
     61    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 
    6362 
    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); 
     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); 
    6867    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    6968 

  • sasmodels/models/triaxial_ellipsoid.py

    r34a9e4e r925ad6e  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
     4All three axes are of different lengths with $R_a \leq R_b \leq R_c$ 
     5**Users should maintain this inequality for all calculations**. 
     6 
     7.. math:: 
     8 
     9    P(q) = \text{scale} V \left< F^2(q) \right> + \text{background} 
     10 
     11where 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 
    418Definition 
    519---------- 
    620 
    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  
    12 Given an ellipsoid 
     21The form factor calculated is 
    1322 
    1423.. math:: 
    1524 
    16     \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1 
    17  
    18 the scattering for randomly oriented particles is defined by the average over 
    19 all orientations $\Omega$ of: 
    20  
    21 .. math:: 
    22  
    23     P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega\Phi^2(qr)\,d\Omega 
    24            + \text{background} 
     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 
    2528 
    2629where 
    … …  
    2831.. math:: 
    2932 
    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 
     33    \Phi(u) = 3 u^{-3} (\sin u - u \cos u) 
    3334 
    34 The $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 
    36 axes, 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  
    45 with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$. 
    46 A 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  
    52 for 
    53  
    54 .. math:: 
    55  
    56     p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1 
    57  
    58 Due 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 
    60 integral 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  
    68 Though for convenience we describe the three radii of the ellipsoid as equatorial 
    69 and polar, they may be given in $any$ size order. To avoid multiple solutions, especially 
    70 with Monte-Carlo fit methods, it may be advisable to restrict their ranges. For typical 
    71 small angle diffraction situations there may be a number of closely similar "best fits", 
    72 so some trial and error, or fixing of some radii at expected values, may help. 
    73      
    7435To provide easy access to the orientation of the triaxial ellipsoid, 
    7536we define the axis of the cylinder using the angles $\theta$, $\phi$ 
    76 and $\psi$. These angles are defined analogously to the elliptical_cylinder below, note that 
    77 angle $\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  
    84 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
    85 see the :ref:`elliptical-cylinder` model for further information. 
     37and $\psi$. These angles are defined on 
     38:numref:`triaxial-ellipsoid-angles` . 
     39The angle $\psi$ is the rotational angle around its own $c$ axis 
     40against the $q$ plane. For example, $\psi = 0$ when the 
     41$a$ axis is parallel to the $x$ axis of the detector. 
    8642 
    8743.. _triaxial-ellipsoid-angles: 
    8844 
    89 .. figure:: img/triaxial_ellipsoid_angle_projection.png 
     45.. figure:: img/triaxial_ellipsoid_angle_projection.jpg 
    9046 
    91     Some examples for an oriented triaxial ellipsoid. 
     47    The angles for oriented ellipsoid. 
    9248 
    9349The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$. 
    9450 
    95 The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the 
     51The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the 
    9652parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major 
    9753equatorial radius, and $R_c$ is the polar radius of the ellipsoid. 
    9854 
    9955NB: The 2nd virial coefficient of the triaxial solid ellipsoid is 
    100 calculated after sorting the three radii to give the most appropriate 
    101 prolate or oblate form, from the new polar radius $R_p = R_c$ and effective equatorial 
    102 radius,  $R_e = \sqrt{R_a R_b}$, to then be used as the effective radius for 
     56calculated based on the polar radius $R_p = R_c$ and equatorial 
     57radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for 
    10358$S(q)$ when $P(q) \cdot S(q)$ is applied. 
    10459 
    … …  
    11469---------- 
    11570 
    116 [1] Finnigan, J.A., Jacobs, D.J., 1971. 
    117 *Light scattering by ellipsoidal particles in solution*, 
    118 J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 
    119  
    120 Authorship 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 
     71L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
     72and Neutron Scattering*, Plenum, New York, 1987. 
    12673""" 
    12774 
    128 from numpy import inf, sin, cos, pi 
     75from numpy import inf 
    12976 
    13077name = "triaxial_ellipsoid" 
    13178title = "Ellipsoid of uniform scattering length density with three independent axes." 
    13279 
    133 description = """ 
    134    Triaxial ellipsoid - see main documentation. 
     80description = """\ 
     81Note: During fitting ensure that the inequality ra<rb<rc is not 
     82        violated. Otherwise the calculation will 
     83        not be correct. 
    13584""" 
    13685category = "shape:ellipsoid" 
    … …  
    14291               "Solvent scattering length density"], 
    14392              ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 
    144                "Minor equatorial radius, Ra"], 
     93               "Minor equatorial radius"], 
    14594              ["radius_equat_major", "Ang", 400, [0, inf], "volume", 
    146                "Major equatorial radius, Rb"], 
     95               "Major equatorial radius"], 
    14796              ["radius_polar", "Ang", 10, [0, inf], "volume", 
    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"], 
     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"], 
    155104             ] 
    156105 
    … …  
    159108def ER(radius_equat_minor, radius_equat_major, radius_polar): 
    160109    """ 
    161     Returns the effective radius used in the S*P calculation 
     110        Returns the effective radius used in the S*P calculation 
    162111    """ 
    163112    import numpy as np 
    164113    from .ellipsoid import ER as ellipsoid_ER 
    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) 
     114    return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major)) 
    174115 
    175116demo = dict(scale=1, background=0, 
    … …  
    183124            phi_pd=15, phi_pd_n=1, 
    184125            psi_pd=15, psi_pd_n=1) 
    185  
    186 q = 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 
    190 qx = q*cos(pi/6.0) 
    191 qy = q*sin(pi/6.0) 
    192 tests = [[{}, 0.05, 24.8839548033], 
    193         [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266 ], 
    194         ] 
    195 del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/two_power_law.py

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

  • sasmodels/rst2html.py

    rf2f5413 r0890871  
    2929from docutils.writers.html4css1 import HTMLTranslator 
    3030from docutils.nodes import SkipNode 
     31 
     32def 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 
     41def 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) 
    3147 
    3248def rst2html(rst, part="whole", math_output="mathjax"): 
    … …  
    139155    assert replace_dollar(u"a (again $in parens$) a") == u"a (again :math:`in parens`) a" 
    140156 
    141 def 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  
    148 def 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  
    157 def view_html_wxapp(html, url=""): 
     157def view_rst_app(filename): 
    158158    import wx  # type: ignore 
    159159    app = wx.App() 
    160     frame = wxview(html, url) 
     160    view_rst(filename) 
    161161    app.MainLoop() 
    162162 
    163 def 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  
    173 def 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  
    185 def 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  
    195 def 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  
    213 def 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) 
    227163 
    228164if __name__ == "__main__": 
    229165    import sys 
    230     view_help(sys.argv[1], qt=True) 
     166    view_rst_app(sys.argv[1]) 
    231167 
Note: See TracChangeset for help on using the changeset viewer.