Changeset f752b9b in sasmodels


Ignore:
Timestamp:
Mar 6, 2019 5:22:31 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket_1156
Children:
a3412a6
Parents:
cc8b183 (diff), e589e9a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'beta_approx' into ticket_1156

Files:
23 edited
6 moved

Legend:

Unmodified
Added
Removed

  • README.rst

    re30d645 r2a64722  
    1010is available. 
    1111 
    12 Example 
     12Install 
    1313------- 
     14 
     15The easiest way to use sasmodels is from `SasView <http://www.sasview.org/>`_. 
     16 
     17You can also install sasmodels as a standalone package in python. Use 
     18`miniconda <https://docs.conda.io/en/latest/miniconda.html>`_ 
     19or `anaconda <https://www.anaconda.com/>`_ 
     20to create a python environment with the sasmodels dependencies:: 
     21 
     22    $ conda create -n sasmodels -c conda-forge numpy scipy matplotlib pyopencl 
     23 
     24The option ``-n sasmodels`` names the environment sasmodels, and the option 
     25``-c conda-forge`` selects the conda-forge package channel because pyopencl 
     26is not part of the base anaconda distribution. 
     27 
     28Activate the environment and install sasmodels:: 
     29 
     30    $ conda activate sasmodels 
     31    (sasmodels) $ pip install sasmodels 
     32 
     33Install `bumps <https://github.com/bumps/bumps>`_ if you want to use it to fit 
     34your data:: 
     35 
     36    (sasmodels) $ pip install bumps 
     37 
     38Usage 
     39----- 
     40 
     41Check that the works:: 
     42 
     43    (sasmodels) $ python -m sasmodels.compare cylinder 
     44 
     45To show the orientation explorer:: 
     46 
     47    (sasmodels) $ python -m sasmodels.jitter 
     48 
     49Documentation is available online as part of the SasView 
     50`fitting perspective <http://www.sasview.org/docs/index.html>`_ 
     51as well as separate pages for 
     52`individual models <http://www.sasview.org/docs/user/sasgui/perspectives/fitting/models/index.html>`_. 
     53Programming details for sasmodels are available in the 
     54`developer documentation <http://www.sasview.org/docs/dev/dev.html>`_. 
     55 
     56 
     57Fitting Example 
     58--------------- 
    1459 
    1560The example directory contains a radial+tangential data set for an oriented 
    1661rod-like shape. 
    1762 
    18 The data is loaded by sas.dataloader from the sasview package, so sasview 
    19 is needed to run the example. 
     63To load the example data, you will need the SAS data loader from the sasview 
     64package. This is not yet available on PyPI, so you will need a copy of the 
     65SasView source code to run it.  Create a directory somewhere to hold the 
     66sasview and sasmodels source code, which we will refer to as $SOURCE. 
    2067 
    21 To run the example, you need sasview, sasmodels and bumps.  Assuming these 
    22 repositories are installed side by side, change to the sasmodels/example 
    23 directory and enter:: 
     68Use the following to install sasview, and the sasmodels examples:: 
    2469 
    25     PYTHONPATH=..:../../sasview/src ../../bumps/run.py fit.py \ 
    26         cylinder --preview 
     70    (sasmodels) $ cd $SOURCE 
     71    (sasmodels) $ conda install git 
     72    (sasmodels) $ git clone https://github.com/sasview/sasview.git 
     73    (sasmodels) $ git clone https://github.com/sasview/sasmodels.git 
    2774 
    28 See bumps documentation for instructions on running the fit.  With the 
    29 python packages installed, e.g., into a virtual environment, then the 
    30 python path need not be set, and the command would be:: 
     75Set the path to the sasview source on your python path within the sasmodels 
     76environment.  On Windows, this will be:: 
    3177 
    32     bumps fit.py cylinder --preview 
     78    (sasmodels)> set PYTHONPATH="$SOURCE\sasview\src" 
     79    (sasmodels)> cd $SOURCE/sasmodels/example 
     80    (sasmodels)> python -m bumps.cli fit.py cylinder --preview 
     81 
     82On Mac/Linux with the standard shell this will be:: 
     83 
     84    (sasmodels) $ export PYTHONPATH="$SOURCE/sasview/src" 
     85    (sasmodels) $ cd $SOURCE/sasmodels/example 
     86    (sasmodels) $ bumps fit.py cylinder --preview 
    3387 
    3488The fit.py model accepts up to two arguments.  The first argument is the 
    … …  
    3892both radial and tangential simultaneously, use the word "both". 
    3993 
    40 Notes 
    41 ----- 
    42  
    43 cylinder.c + cylinder.py is the cylinder model with renamed variables and 
    44 sld scaled by 1e6 so the numbers are nicer.  The model name is "cylinder" 
    45  
    46 lamellar.py is an example of a single file model with embedded C code. 
     94See `bumps documentation <https://bumps.readthedocs.io/>`_ for detailed 
     95instructions on running the fit. 
    4796 
    4897|TravisStatus|_ 

  • explore/beta/sasfit_compare.py

    r2a12351b r119073a  
    505505    } 
    506506 
    507     Q, IQ = load_sasfit(data_file('richard_test.txt')) 
    508     Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
    509     Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     507    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt')) 
     508    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt')) 
     509    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt')) 
    510510    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    511511    actual = sphere_r(Q, norm="sasfit", **pars) 
    … …  
    526526    } 
    527527 
    528     Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
    529     Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
    530     Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     528    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQD.txt')) 
     529    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQSD.txt')) 
     530    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQBD.txt')) 
    531531    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    532532    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 

  • explore/precision.py

    raa8c6e0 rcd28947  
    207207    return model_info 
    208208 
    209 # Hack to allow second parameter A in two parameter functions 
     209# Hack to allow second parameter A in the gammainc and gammaincc functions. 
     210# Create a 2-D variant of the precision test if we need to handle other two 
     211# parameter functions. 
    210212A = 1 
    211213def parse_extra_pars(): 
     214    """ 
     215    Parse the command line looking for the second parameter "A=..." for the 
     216    gammainc/gammaincc functions. 
     217    """ 
    212218    global A 
    213219 
    … …  
    333339) 
    334340add_function( 
     341    # Note: "a" is given as A=... on the command line via parse_extra_pars 
    335342    name="gammainc(x)", 
    336343    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
    … …  
    339346) 
    340347add_function( 
     348    # Note: "a" is given as A=... on the command line via parse_extra_pars 
    341349    name="gammaincc(x)", 
    342350    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 

  • sasmodels/compare.py

    r07646b6 rc1799d3  
    11521152        'rel_err'   : True, 
    11531153        'explore'   : False, 
    1154         'use_demo'  : True, 
     1154        'use_demo'  : False, 
    11551155        'zero'      : False, 
    11561156        'html'      : False, 

  • sasmodels/direct_model.py

    r5024a56 rc1799d3  
    332332 
    333333        # Need to pull background out of resolution for multiple scattering 
    334         background = pars.get('background', DEFAULT_BACKGROUND) 
     334        default_background = self._model.info.parameters.common_parameters[1].default 
     335        background = pars.get('background', default_background) 
    335336        pars = pars.copy() 
    336337        pars['background'] = 0. 

  • sasmodels/generate.py

    r39a06c9 ra8a1d48  
    703703    """ 
    704704    for code in source: 
    705         m = _FQ_PATTERN.search(code) 
    706         if m is not None: 
     705        if _FQ_PATTERN.search(code) is not None: 
    707706            return True 
    708707    return False 
    … …  
    712711    # type: (List[str]) -> bool 
    713712    """ 
    714     Return True if C source defines "void Fq(". 
     713    Return True if C source defines "double shell_volume(". 
    715714    """ 
    716715    for code in source: 
    717         m = _SHELL_VOLUME_PATTERN.search(code) 
    718         if m is not None: 
     716        if _SHELL_VOLUME_PATTERN.search(code) is not None: 
    719717            return True 
    720718    return False 
    … …  
    10081006        pars = model_info.parameters.kernel_parameters 
    10091007    else: 
    1010         pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     1008        pars = (model_info.parameters.common_parameters 
     1009                + model_info.parameters.kernel_parameters) 
    10111010    partable = make_partable(pars) 
    10121011    subst = dict(id=model_info.id.replace('_', '-'), 

  • sasmodels/jitter.py

    r1198f90 r7d97437  
    1515    pass 
    1616 
     17import matplotlib as mpl 
    1718import matplotlib.pyplot as plt 
    1819from matplotlib.widgets import Slider 
    … …  
    746747        pass 
    747748 
    748     axcolor = 'lightgoldenrodyellow' 
     749    # CRUFT: use axisbg instead of facecolor for matplotlib<2 
     750    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg' 
     751    props = {facecolor_prop: 'lightgoldenrodyellow'} 
    749752 
    750753    ## add control widgets to plot 
    751     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    752     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    753     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
     754    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) 
     755    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) 
     756    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) 
    754757    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    755758    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    756759    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    757760 
    758     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    759     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    760     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     761    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) 
     762    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) 
     763    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) 
    761764    # Note: using ridiculous definition of rectangle distribution, whose width 
    762765    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 

  • sasmodels/kernel.py

    re44432d rcd28947  
    133133        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
    134134        total_weight = self.result[nout*self.q_input.nq + 0] 
     135        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it 
     136        # is okay to test directly against zero.  If weight is zero then I(q), 
     137        # etc. must also be zero. 
    135138        if total_weight == 0.: 
    136139            total_weight = 1. 

  • sasmodels/kernelcl.py

    rf872fd1 r8064d5e  
    5858import time 
    5959 
     60try: 
     61    from time import perf_counter as clock 
     62except ImportError: # CRUFT: python < 3.3 
     63    import sys 
     64    if sys.platform.count("darwin") > 0: 
     65        from time import time as clock 
     66    else: 
     67        from time import clock 
     68 
    6069import numpy as np  # type: ignore 
    61  
    6270 
    6371# Attempt to setup opencl. This may fail if the pyopencl package is not 
    … …  
    611619        #Call kernel and retrieve results 
    612620        wait_for = None 
    613         last_nap = time.clock() 
     621        last_nap = clock() 
    614622        step = 1000000//self.q_input.nq + 1 
    615623        for start in range(0, call_details.num_eval, step): 
    … …  
    622630                # Allow other processes to run 
    623631                wait_for[0].wait() 
    624                 current_time = time.clock() 
     632                current_time = clock() 
    625633                if current_time - last_nap > 0.5: 
    626634                    time.sleep(0.001) 

  • sasmodels/modelinfo.py

    r39a06c9 rc1799d3  
    404404      parameters counted as n individual parameters p1, p2, ... 
    405405 
     406    * *common_parameters* is the list of common parameters, with a unique 
     407      copy for each model so that structure factors can have a default 
     408      background of 0.0. 
     409 
    406410    * *call_parameters* is the complete list of parameters to the kernel, 
    407411      including scale and background, with vector parameters recorded as 
    … …  
    416420    parameters don't use vector notation, and instead use p1, p2, ... 
    417421    """ 
    418     # scale and background are implicit parameters 
    419     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    420  
    421422    def __init__(self, parameters): 
    422423        # type: (List[Parameter]) -> None 
     424 
     425        # scale and background are implicit parameters 
     426        # Need them to be unique to each model in case they have different 
     427        # properties, such as default=0.0 for structure factor backgrounds. 
     428        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] 
     429 
    423430        self.kernel_parameters = parameters 
    424431        self._set_vector_lengths() 
    … …  
    468475                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    469476        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     477 
     478    def set_zero_background(self): 
     479        """ 
     480        Set the default background to zero for this model.  This is done for 
     481        structure factor models. 
     482        """ 
     483        # type: () -> None 
     484        # Make sure background is the second common parameter. 
     485        assert self.common_parameters[1].id == "background" 
     486        self.common_parameters[1].default = 0.0 
     487        self.defaults = self._get_defaults() 
    470488 
    471489    def check_angles(self): 
    … …  
    567585    def _get_call_parameters(self): 
    568586        # type: () -> List[Parameter] 
    569         full_list = self.COMMON[:] 
     587        full_list = self.common_parameters[:] 
    570588        for p in self.kernel_parameters: 
    571589            if p.length == 1: 
    … …  
    670688 
    671689        # Gather the user parameters in order 
    672         result = control + self.COMMON 
     690        result = control + self.common_parameters 
    673691        for p in self.kernel_parameters: 
    674692            if not is2d and p.type in ('orientation', 'magnetic'): 
    … …  
    770788 
    771789    info = ModelInfo() 
     790 
     791    # Build the parameter table 
    772792    #print("make parameter table", kernel_module.parameters) 
    773793    parameters = make_parameter_table(getattr(kernel_module, 'parameters', [])) 
     794 
     795    # background defaults to zero for structure factor models 
     796    structure_factor = getattr(kernel_module, 'structure_factor', False) 
     797    if structure_factor: 
     798        parameters.set_zero_background() 
     799 
     800    # TODO: remove demo parameters 
     801    # The plots in the docs are generated from the model default values. 
     802    # Sascomp set parameters from the command line, and so doesn't need 
     803    # demo values for testing. 
    774804    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None)) 
     805 
    775806    filename = abspath(kernel_module.__file__).replace('.pyc', '.py') 
    776807    kernel_id = splitext(basename(filename))[0] 

  • sasmodels/models/hardsphere.py

    r304c775 rc1799d3  
    162162    return pars 
    163163 
    164 demo = dict(radius_effective=200, volfraction=0.2, 
    165             radius_effective_pd=0.1, radius_effective_pd_n=40) 
    166164# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1, 
    167165# assuming double precision sasview is correct 

  • sasmodels/models/pearl_necklace.c

    r99658f6 r9b5fd42  
    4040    const double si = sas_sinx_x(q*A_s); 
    4141    const double omsi = 1.0 - si; 
    42     const double pow_si = pow(si, num_pearls); 
     42    const double pow_si = pown(si, num_pearls); 
    4343 
    4444    // form factor for num_pearls 
    … …  
    8181radius_from_volume(double radius, double edge_sep, double thick_string, double fp_num_pearls) 
    8282{ 
    83     const int num_pearls = (int) fp_num_pearls +0.5; 
    8483    const double vol_tot = form_volume(radius, edge_sep, thick_string, fp_num_pearls); 
    8584    return cbrt(vol_tot/M_4PI_3); 

  • sasmodels/resolution.py

    r9e7837a re2592f0  
    445445    q = np.sort(q) 
    446446    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
    447         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
     447        n_low = int(np.ceil((q[0]-q_min) / (q[1]-q[0]))) if q[1] > q[0] else 15 
    448448        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    449449    else: 
    450450        q_low = [] 
    451451    if q_max - 2*MINIMUM_RESOLUTION > q[-1]: 
    452         n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 
     452        n_high = int(np.ceil((q_max-q[-1]) / (q[-1]-q[-2]))) if q[-1] > q[-2] else 15 
    453453        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    454454    else: 
    … …  
    499499            q_min = q[0]*MINIMUM_ABSOLUTE_Q 
    500500        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    501         q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
     501        q_low = np.logspace(log10(q_min), log10(q[0]), int(np.ceil(n_low))+1)[:-1] 
    502502    else: 
    503503        q_low = [] 
    504504    if q_max > q[-1]: 
    505505        n_high = log_delta_q * (log(q_max)-log(q[-1])) 
    506         q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:] 
     506        q_high = np.logspace(log10(q[-1]), log10(q_max), int(np.ceil(n_high))+1)[1:] 
    507507    else: 
    508508        q_high = [] 

  • sasmodels/sasview_model.py

    r5024a56 ra8a1d48  
    382382            hidden.add('scale') 
    383383            hidden.add('background') 
    384             self._model_info.parameters.defaults['background'] = 0. 
    385384 
    386385        # Update the parameter lists to exclude any hidden parameters 
    … …  
    695694            return self._calculate_Iq(qx, qy) 
    696695 
    697     def _calculate_Iq(self, qx, qy=None, Fq=False, effective_radius_type=1): 
     696    def _calculate_Iq(self, qx, qy=None): 
    698697        if self._model is None: 
    699698            self._model = core.build_model(self._model_info) 
    … …  
    715714        #print("values", values) 
    716715        #print("is_mag", is_magnetic) 
    717         if Fq: 
    718             result = calculator.Fq(call_details, values, cutoff=self.cutoff, 
    719                                    magnetic=is_magnetic, 
    720                                    effective_radius_type=effective_radius_type) 
    721716        result = calculator(call_details, values, cutoff=self.cutoff, 
    722717                            magnetic=is_magnetic) 
    … …  
    736731        Calculate the effective radius for P(q)*S(q) 
    737732 
     733        *mode* is the R_eff type, which defaults to 1 to match the ER 
     734        calculation for sasview models from version 3.x. 
     735 
    738736        :return: the value of the effective radius 
    739737        """ 
    740         Fq = self._calculate_Iq([0.1], True, mode) 
    741         return Fq[2] 
     738        # ER and VR are only needed for old multiplication models, based on 
     739        # sas.sascalc.fit.MultiplicationModel.  Fail for now.  If we want to 
     740        # continue supporting them then add some test cases so that the code 
     741        # is exercised.  We can access ER/VR using the kernel Fq function by 
     742        # extending _calculate_Iq so that it calls: 
     743        #    if er_mode > 0: 
     744        #        res = calculator.Fq(call_details, values, cutoff=self.cutoff, 
     745        #                            magnetic=False, effective_radius_type=mode) 
     746        #        R_eff, form_shell_ratio = res[2], res[4] 
     747        #        return R_eff, form_shell_ratio 
     748        # Then use the following in calculate_ER: 
     749        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=mode) 
     750        #    return ER 
     751        # Similarly, for calculate_VR: 
     752        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=1) 
     753        #    return VR 
     754        # Obviously a combined calculate_ER_VR method would be better, but 
     755        # we only need them to support very old models, so ignore the 2x 
     756        # performance hit. 
     757        raise NotImplementedError("ER function is no longer available.") 
    742758 
    743759    def calculate_VR(self): 
    … …  
    748764        :return: the value of the form:shell volume ratio 
    749765        """ 
    750         Fq = self._calculate_Iq([0.1], True, mode) 
    751         return Fq[4] 
     766        # See comments in calculate_ER. 
     767        raise NotImplementedError("VR function is no longer available.") 
    752768 
    753769    def set_dispersion(self, parameter, dispersion): 
    … …  
    914930    CylinderModel().evalDistribution([0.1, 0.1]) 
    915931 
     932def test_structure_factor_background(): 
     933    # type: () -> None 
     934    """ 
     935    Check that sasview model and direct model match, with background=0. 
     936    """ 
     937    from .data import empty_data1D 
     938    from .core import load_model_info, build_model 
     939    from .direct_model import DirectModel 
     940 
     941    model_name = "hardsphere" 
     942    q = [0.0] 
     943 
     944    sasview_model = _make_standard_model(model_name)() 
     945    sasview_value = sasview_model.evalDistribution(np.array(q))[0] 
     946 
     947    data = empty_data1D(q) 
     948    model_info = load_model_info(model_name) 
     949    model = build_model(model_info) 
     950    direct_model = DirectModel(data, model) 
     951    direct_value_zero_background = direct_model(background=0.0) 
     952 
     953    assert sasview_value == direct_value_zero_background 
     954 
     955    # Additionally check that direct value background defaults to zero 
     956    direct_value_default = direct_model() 
     957    assert sasview_value == direct_value_default 
     958 
     959 
    916960def magnetic_demo(): 
    917961    Model = _make_standard_model('sphere') 
    … …  
    934978    #print("rpa:", test_rpa()) 
    935979    #test_empty_distribution() 
     980    #test_structure_factor_background() 

  • sasmodels/weights.py

    r3d58247 re2592f0  
    2323    default = dict(npts=35, width=0, nsigmas=3) 
    2424    def __init__(self, npts=None, width=None, nsigmas=None): 
    25         self.npts = self.default['npts'] if npts is None else npts 
     25        self.npts = self.default['npts'] if npts is None else int(npts) 
    2626        self.width = self.default['width'] if width is None else width 
    2727        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas 

  • explore/realspace.py

    r362a66f r5778279  
    99import numpy as np 
    1010from numpy import pi, radians, sin, cos, sqrt 
    11 from numpy.random import poisson, uniform, randn, rand 
     11from numpy.random import poisson, uniform, randn, rand, randint 
    1212from numpy.polynomial.legendre import leggauss 
    1313from scipy.integrate import simps 
    … …  
    7878 
    7979 
     80I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 
     81 
    8082class Shape: 
    81     rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 
     83    rotation = I3 
    8284    center = np.array([0., 0., 0.])[:, None] 
    8385    r_max = None 
     86    lattice_size = np.array((1, 1, 1)) 
     87    lattice_spacing = np.array((1., 1., 1.)) 
     88    lattice_distortion = 0.0 
     89    lattice_rotation = 0.0 
     90    lattice_type = "" 
    8491 
    8592    def volume(self): 
    … …  
    96103 
    97104    def rotate(self, theta, phi, psi): 
    98         self.rotation = rotation(theta, phi, psi) * self.rotation 
     105        if theta != 0. or phi != 0. or psi != 0.: 
     106            self.rotation = rotation(theta, phi, psi) * self.rotation 
    99107        return self 
    100108 
    … …  
    103111        return self 
    104112 
     113    def lattice(self, size=(1, 1, 1), spacing=(2, 2, 2), type="sc", 
     114                distortion=0.0, rotation=0.0): 
     115        self.lattice_size = np.asarray(size, 'i') 
     116        self.lattice_spacing = np.asarray(spacing, 'd') 
     117        self.lattice_type = type 
     118        self.lattice_distortion = distortion 
     119        self.lattice_rotation = rotation 
     120 
    105121    def _adjust(self, points): 
    106         points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 
     122        if self.rotation is I3: 
     123            points = points.T + self.center 
     124        else: 
     125            points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 
     126        if self.lattice_type: 
     127            points = self._apply_lattice(points) 
    107128        return points.T 
    108129 
    109     def r_bins(self, q, over_sampling=1, r_step=0.): 
    110         r_max = min(2 * pi / q[0], self.r_max) 
     130    def r_bins(self, q, over_sampling=10, r_step=0.): 
     131        if self.lattice_type: 
     132            r_max = np.sqrt(np.sum(self.lattice_size*self.lattice_spacing*self.dims)**2)/2 
     133        else: 
     134            r_max = self.r_max 
     135        #r_max = min(2 * pi / q[0], r_max) 
    111136        if r_step == 0.: 
    112137            r_step = 2 * pi / q[-1] / over_sampling 
    113138        #r_step = 0.01 
    114139        return np.arange(r_step, r_max, r_step) 
     140 
     141    def _apply_lattice(self, points): 
     142        """Spread points to different lattice positions""" 
     143        size = self.lattice_size 
     144        spacing = self.lattice_spacing 
     145        shuffle = self.lattice_distortion 
     146        rotate = self.lattice_rotation 
     147        lattice = self.lattice_type 
     148 
     149        if rotate != 0: 
     150            # To vectorize the rotations we will need to unwrap the matrix multiply 
     151            raise NotImplementedError("don't handle rotations yet") 
     152 
     153        # Determine the number of lattice points in the lattice 
     154        shapes_per_cell = 2 if lattice == "bcc" else 4 if lattice == "fcc" else 1 
     155        number_of_lattice_points = np.prod(size) * shapes_per_cell 
     156 
     157        # For each point in the original shape, figure out which lattice point 
     158        # to translate it to.  This is both cell index (i*ny*nz + j*nz  + k) as 
     159        # well as the point in the cell (corner, body center or face center). 
     160        nsamples = points.shape[1] 
     161        lattice_point = randint(number_of_lattice_points, size=nsamples) 
     162 
     163        # Translate the cell index into the i,j,k coordinates of the senter 
     164        cell_index = lattice_point // shapes_per_cell 
     165        center = np.vstack((cell_index//(size[1]*size[2]), 
     166                            (cell_index%(size[1]*size[2]))//size[2], 
     167                            cell_index%size[2])) 
     168        center = np.asarray(center, dtype='d') 
     169        if lattice == "bcc": 
     170            center[:, lattice_point % shapes_per_cell == 1] += [[0.5], [0.5], [0.5]] 
     171        elif lattice == "fcc": 
     172            center[:, lattice_point % shapes_per_cell == 1] += [[0.0], [0.5], [0.5]] 
     173            center[:, lattice_point % shapes_per_cell == 2] += [[0.5], [0.0], [0.5]] 
     174            center[:, lattice_point % shapes_per_cell == 3] += [[0.5], [0.5], [0.0]] 
     175 
     176        # Each lattice point has its own displacement from the ideal position. 
     177        # Not checking that shapes do not overlap if displacement is too large. 
     178        offset = shuffle*(randn(3, number_of_lattice_points) if shuffle < 0.3 
     179                          else rand(3, number_of_lattice_points)) 
     180        center += offset[:, cell_index] 
     181 
     182        # Each lattice point has its own rotation.  Rotate the point prior to 
     183        # applying any displacement. 
     184        # rotation = rotate*(randn(size=(shapes, 3)) if shuffle < 30 else rand(size=(nsamples, 3))) 
     185        # for k in shapes: points[k] = rotation[k]*points[k] 
     186        points += center*(np.array([spacing])*np.array(self.dims)).T 
     187        return points 
    115188 
    116189class Composite(Shape): 
    … …  
    669742    Iq = 100 * np.ones_like(qx) 
    670743    data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 
    671     data.x_bins = qx[0,:] 
    672     data.y_bins = qy[:,0] 
     744    data.x_bins = qx[0, :] 
     745    data.y_bins = qy[:, 0] 
    673746    data.filename = "fake data" 
    674747 
    … …  
    695768    return shape, fn, fn_xy 
    696769 
    697 def build_sphere(radius=125, rho=2): 
     770DEFAULT_SPHERE_RADIUS = 125 
     771DEFAULT_SPHERE_CONTRAST = 2 
     772def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST): 
    698773    shape = TriaxialEllipsoid(radius, radius, radius, rho) 
    699774    fn = lambda q: sphere_Iq(q, radius)*rho**2 
    … …  
    751826    return shape, fn, fn_xy 
    752827 
    753 def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
    754                   shuffle=0, rotate=0): 
     828def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     829                        shuffle=0, rotate=0): 
    755830    a, b, c = shape.dims 
    756     shapes = [copy(shape) 
     831    corners= [copy(shape) 
    757832              .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
    758833                     (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
    … …  
    762837              for iy in range(ny) 
    763838              for iz in range(nz)] 
    764     lattice = Composite(shapes) 
     839    lattice = Composite(corners) 
    765840    return lattice 
    766841 
     842def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     843                      shuffle=0, rotate=0): 
     844    a, b, c = shape.dims 
     845    corners = [copy(shape) 
     846               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     847                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     848                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     849               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     850               for ix in range(nx) 
     851               for iy in range(ny) 
     852               for iz in range(nz)] 
     853    centers = [copy(shape) 
     854               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     855                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     856                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     857               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     858               for ix in range(nx) 
     859               for iy in range(ny) 
     860               for iz in range(nz)] 
     861    lattice = Composite(corners + centers) 
     862    return lattice 
     863 
     864def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     865                      shuffle=0, rotate=0): 
     866    a, b, c = shape.dims 
     867    corners = [copy(shape) 
     868               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     869                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     870                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     871               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     872               for ix in range(nx) 
     873               for iy in range(ny) 
     874               for iz in range(nz)] 
     875    faces_a = [copy(shape) 
     876               .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     877                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     878                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     879               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     880               for ix in range(nx) 
     881               for iy in range(ny) 
     882               for iz in range(nz)] 
     883    faces_b = [copy(shape) 
     884               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     885                      (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     886                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     887               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     888               for ix in range(nx) 
     889               for iy in range(ny) 
     890               for iz in range(nz)] 
     891    faces_c = [copy(shape) 
     892               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     893                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     894                      (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     895               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     896               for ix in range(nx) 
     897               for iy in range(ny) 
     898               for iz in range(nz)] 
     899    lattice = Composite(corners + faces_a + faces_b + faces_c) 
     900    return lattice 
    767901 
    768902SHAPE_FUNCTIONS = OrderedDict([ 
    … …  
    775909]) 
    776910SHAPES = list(SHAPE_FUNCTIONS.keys()) 
     911LATTICE_FUNCTIONS = OrderedDict([ 
     912    ("sc", build_sc_lattice), 
     913    ("bcc", build_bcc_lattice), 
     914    ("fcc", build_fcc_lattice), 
     915]) 
     916LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys()) 
    777917 
    778918def check_shape(title, shape, fn=None, show_points=False, 
    … …  
    783923    r = shape.r_bins(q, r_step=r_step) 
    784924    sampling_density = samples / shape.volume 
     925    print("sampling points") 
    785926    rho, points = shape.sample(sampling_density) 
     927    print("calculating Pr") 
    786928    t0 = time.time() 
    787929    Pr = calc_Pr(r, rho-rho_solvent, points) 
    … …  
    792934    import pylab 
    793935    if show_points: 
    794          plot_points(rho, points); pylab.figure() 
     936        plot_points(rho, points); pylab.figure() 
    795937    plot_calc(r, Pr, q, Iq, theory=theory, title=title) 
    796938    pylab.gcf().canvas.set_window_title(title) 
    … …  
    806948    Qx, Qy = np.meshgrid(qx, qy) 
    807949    sampling_density = samples / shape.volume 
     950    print("sampling points") 
    808951    t0 = time.time() 
    809952    rho, points = shape.sample(sampling_density) 
    … …  
    844987                        help='lattice size') 
    845988    parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 
    846                         help='lattice spacing') 
     989                        help='lattice spacing (relative to shape)') 
     990    parser.add_argument('-t', '--type', choices=LATTICE_TYPES, 
     991                        default=LATTICE_TYPES[0], 
     992                        help='lattice type') 
    847993    parser.add_argument('-r', '--rotate', type=float, default=0., 
    848994                        help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 
    … …  
    8581004    nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 
    8591005    dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 
    860     shuffle, rotate = opts.shuffle, opts.rotate 
     1006    distortion, rotation = opts.shuffle, opts.rotate 
    8611007    shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 
    862     if nx > 1 or ny > 1 or nz > 1: 
    863         shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     1008    view = tuple(float(v) for v in opts.view.split(',')) 
     1009    # If comparing a sphere in a cubic lattice, compare against the 
     1010    # corresponding paracrystalline model. 
     1011    if opts.shape == "sphere" and dx == dy == dz and nx*ny*nz > 1: 
     1012        radius = pars.get('radius', DEFAULT_SPHERE_RADIUS) 
     1013        model_name = opts.type + "_paracrystal" 
     1014        model_pars = { 
     1015            "scale": 1., 
     1016            "background": 0., 
     1017            "lattice_spacing": 2*radius*dx, 
     1018            "lattice_distortion": distortion, 
     1019            "radius": radius, 
     1020            "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST), 
     1021            "sld_solvent": 0., 
     1022            "theta": view[0], 
     1023            "phi": view[1], 
     1024            "psi": view[2], 
     1025        } 
     1026        fn, fn_xy = wrap_sasmodel(model_name, **model_pars) 
     1027    if nx*ny*nz > 1: 
     1028        if rotation != 0: 
     1029            print("building %s lattice"%opts.type) 
     1030            build_lattice = LATTICE_FUNCTIONS[opts.type] 
     1031            shape = build_lattice(shape, nx, ny, nz, dx, dy, dz, 
     1032                                  distortion, rotation) 
     1033        else: 
     1034            shape.lattice(size=(nx, ny, nz), spacing=(dx, dy, dz), 
     1035                          type=opts.type, 
     1036                          rotation=rotation, distortion=distortion) 
     1037 
    8641038    title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 
    8651039    if opts.dim == 1: 
    … …  
    8671041                    mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
    8681042    else: 
    869         view = tuple(float(v) for v in opts.view.split(',')) 
    8701043        check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 
    8711044                       mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 

  • sasmodels/convert.py

    r610ef23 rb3f4831  
    166166        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 
    167167    if version < (4, 2, 0): 
     168        oldpars = _hand_convert_4_1_to_4_2(name, oldpars) 
    168169        oldpars = _rename_magnetic_pars(oldpars) 
     170    return oldpars 
     171 
     172def _hand_convert_4_1_to_4_2(name, oldpars): 
     173    if name in ('bcc_paracrystal', 'fcc_paracrystal', 'sc_paracrystal'): 
     174        oldpars['lattice_spacing'] = oldpars.pop('dnn') 
     175        oldpars['lattice_distortion'] = oldpars.pop('d_factor') 
    169176    return oldpars 
    170177 

  • sasmodels/models/bcc_paracrystal.c

    r108e70e r642046e  
    11static double 
    2 bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     2bcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 
    33{ 
    44    // Equations from Matsuoka 26-27-28, multiplied by |q| 
    … …  
    1717    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
    1818    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
    19     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     19    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    2020    const double exp_arg = exp(arg); 
    2121    const double Zq = -cube(expm1(2.0*arg)) 
    22         / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
    23           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
    24           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     22        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 
    2525 
    2626#elif 0 
    … …  
    3636    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)] 
    3737    // 
    38     const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     38    const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    3939    const double sinh_qd = sinh(arg); 
    4040    const double cosh_qd = cosh(arg); 
    41     const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
    42                     * sinh_qd/(cosh_qd - cos(dnn*a2)) 
    43                     * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     41    const double Zq = sinh_qd/(cosh_qd - cos(lattice_spacing*a1)) 
     42                    * sinh_qd/(cosh_qd - cos(lattice_spacing*a2)) 
     43                    * sinh_qd/(cosh_qd - cos(lattice_spacing*a3)); 
    4444#else 
    45     const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     45    const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    4646    const double tanh_qd = tanh(arg); 
    4747    const double cosh_qd = cosh(arg); 
    48     const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd) 
    49                     * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd) 
    50                     * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd); 
     48    const double Zq = tanh_qd/(1.0 - cos(lattice_spacing*a1)/cosh_qd) 
     49                    * tanh_qd/(1.0 - cos(lattice_spacing*a2)/cosh_qd) 
     50                    * tanh_qd/(1.0 - cos(lattice_spacing*a3)/cosh_qd); 
    5151#endif 
    5252 
    … …  
    5757// occupied volume fraction calculated from lattice symmetry and sphere radius 
    5858static double 
    59 bcc_volume_fraction(double radius, double dnn) 
     59bcc_volume_fraction(double radius, double lattice_spacing) 
    6060{ 
    61     return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
     61    return 2.0*sphere_volume(radius/lattice_spacing); 
    6262} 
    6363 
    … …  
    6969 
    7070 
    71 static double Iq(double q, double dnn, 
    72     double d_factor, double radius, 
     71static double Iq(double q, double lattice_spacing, 
     72    double lattice_distortion, double radius, 
    7373    double sld, double solvent_sld) 
    7474{ 
    … …  
    9494            const double qa = qab*cos_phi; 
    9595            const double qb = qab*sin_phi; 
    96             const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     96            const double form = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    9797            inner_sum += GAUSS_W[j] * form; 
    9898        } 
    … …  
    103103    const double Zq = outer_sum/(4.0*M_PI); 
    104104    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    105     return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
     105    return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    106106} 
    107107 
    108108 
    109109static double Iqabc(double qa, double qb, double qc, 
    110     double dnn, double d_factor, double radius, 
     110    double lattice_spacing, double lattice_distortion, double radius, 
    111111    double sld, double solvent_sld) 
    112112{ 
    113113    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
    114     const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     114    const double Zq = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    115115    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    116     return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
     116    return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    117117} 

  • sasmodels/models/bcc_paracrystal.py

    rda7b26b rb3f4831  
    3434.. math:: 
    3535 
    36     V_\text{lattice} = \frac{16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3} 
     36    V_\text{lattice} = \frac{8\pi}{3} \frac{R^3}{\left(2D/\sqrt{3}\right)^3} 
    3737 
    3838 
    … …  
    104104 
    105105* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    106 * **Last Modified by:** Paul Butler **Date:** September 29, 2016 
    107 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 
     106* **Last Modified by:** Paul Butler **Date:** September 16, 2018 
     107* **Last Reviewed by:** Paul Butler **Date:** September 16, 2018 
    108108""" 
    109109 
    … …  
    127127# pylint: disable=bad-whitespace, line-too-long 
    128128#             ["name", "units", default, [lower, upper], "type","description" ], 
    129 parameters = [["dnn",         "Ang",       220,    [-inf, inf], "",            "Nearest neighbour distance"], 
    130               ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"], 
     129parameters = [["lattice_spacing",         "Ang",       220,    [-inf, inf], "",            "Lattice spacing"], 
     130              ["lattice_distortion",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"], 
    131131              ["radius",      "Ang",        40,    [0, inf],    "volume",      "Particle radius"], 
    132132              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "sld",         "Particle scattering length density"], 
    … …  
    149149    # in this range 'cuz its easy. 
    150150    radius = 10**np.random.uniform(1.3, 4) 
    151     d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
    152     dnn_fraction = np.random.beta(a=10, b=1) 
    153     dnn = radius*4/np.sqrt(3)/dnn_fraction 
     151    lattice_distortion = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
     152    lattice_spacing_fraction = np.random.beta(a=10, b=1) 
     153    lattice_spacing = radius*4/np.sqrt(3)/lattice_spacing_fraction 
    154154    pars = dict( 
    155155        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
    156         dnn=dnn, 
    157         d_factor=d_factor, 
     156        lattice_spacing=lattice_spacing, 
     157        lattice_distortion=lattice_distortion, 
    158158        radius=radius, 
    159159    ) 

  • sasmodels/models/fcc_paracrystal.c

    r71b751d rcc8b183  
    11static double 
    2 fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     2fcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 
    33{ 
    44    // Equations from Matsuoka 17-18-19, multiplied by |q| 
    … …  
    1616    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
    1717    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
    18     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     18    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    1919    const double exp_arg = exp(arg); 
    2020    const double Zq = -cube(expm1(2.0*arg)) 
    21         / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
    22           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
    23           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     21        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 
    2424 
    2525    return Zq; 
    … …  
    2929// occupied volume fraction calculated from lattice symmetry and sphere radius 
    3030static double 
    31 fcc_volume_fraction(double radius, double dnn) 
     31fcc_volume_fraction(double radius, double lattice_spacing) 
    3232{ 
    33     return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
     33    return 4.0*sphere_volume(radius/lattice_spacing); 
    3434} 
    3535 
    … …  
    4141 
    4242 
    43 static double Iq(double q, double dnn, 
    44   double d_factor, double radius, 
     43static double Iq(double q, double lattice_spacing, 
     44  double lattice_distortion, double radius, 
    4545  double sld, double solvent_sld) 
    4646{ 
    … …  
    6666            const double qa = qab*cos_phi; 
    6767            const double qb = qab*sin_phi; 
    68             const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     68            const double form = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    6969            inner_sum += GAUSS_W[j] * form; 
    7070        } 
    … …  
    7676    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    7777 
    78     return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     78    return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    7979} 
    8080 
    8181static double Iqabc(double qa, double qb, double qc, 
    82     double dnn, double d_factor, double radius, 
     82    double lattice_spacing, double lattice_distortion, double radius, 
    8383    double sld, double solvent_sld) 
    8484{ 
    8585    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
    8686    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    87     const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
    88     return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     87    const double Zq = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
     88    return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    8989} 

  • sasmodels/models/fcc_paracrystal.py

    rda7b26b rb3f4831  
    33#note - calculation requires double precision 
    44r""" 
    5 .. warning:: This model and this model description are under review following  
    6              concerns raised by SasView users. If you need to use this model,  
    7              please email help@sasview.org for the latest situation. *The  
     5.. warning:: This model and this model description are under review following 
     6             concerns raised by SasView users. If you need to use this model, 
     7             please email help@sasview.org for the latest situation. *The 
    88             SasView Developers. September 2018.* 
    99 
    … …  
    100100 
    101101Authorship and Verification 
    102 --------------------------- 
     102---------------------------- 
    103103 
    104104* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    105 * **Last Modified by:** Paul Butler **Date:** September 29, 2016 
    106 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 
     105* **Last Modified by:** Paul Butler **Date:** September 16, 2018 
     106* **Last Reviewed by:** Paul Butler **Date:** September 16, 2018 
    107107""" 
    108108 
    … …  
    123123# pylint: disable=bad-whitespace, line-too-long 
    124124#             ["name", "units", default, [lower, upper], "type","description"], 
    125 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], 
    126               ["d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 
     125parameters = [["lattice_spacing", "Ang", 220, [-inf, inf], "", "Lattice spacing"], 
     126              ["lattice_distortion", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 
    127127              ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"], 
    128128              ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], 
    … …  
    139139    # copied from bcc_paracrystal 
    140140    radius = 10**np.random.uniform(1.3, 4) 
    141     d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
    142     dnn_fraction = np.random.beta(a=10, b=1) 
    143     dnn = radius*4/np.sqrt(2)/dnn_fraction 
     141    lattice_distortion = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
     142    lattice_spacing_fraction = np.random.beta(a=10, b=1) 
     143    lattice_spacing = radius*4/np.sqrt(2)/lattice_spacing_fraction 
    144144    pars = dict( 
    145145        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
    146         dnn=dnn, 
    147         d_factor=d_factor, 
     146        latice_spacing=lattice_spacing, 
     147        lattice_distortion=d_factor, 
    148148        radius=radius, 
    149149    ) 

  • sasmodels/models/sc_paracrystal.c

    r71b751d rcc8b183  
    11static double 
    2 sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     2sc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 
    33{ 
    44    // Equations from Matsuoka 9-10-11, multiplied by |q| 
    … …  
    1616    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
    1717    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
    18     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     18    const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 
    1919    const double exp_arg = exp(arg); 
    2020    const double Zq = -cube(expm1(2.0*arg)) 
    21         / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
    22           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
    23           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     21        / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 
    2424 
    2525    return Zq; 
    … …  
    2828// occupied volume fraction calculated from lattice symmetry and sphere radius 
    2929static double 
    30 sc_volume_fraction(double radius, double dnn) 
     30sc_volume_fraction(double radius, double lattice_spacing) 
    3131{ 
    32     return sphere_volume(radius/dnn); 
     32    return sphere_volume(radius/lattice_spacing); 
    3333} 
    3434 
    … …  
    4141 
    4242static double 
    43 Iq(double q, double dnn, 
    44     double d_factor, double radius, 
     43Iq(double q, double lattice_spacing, 
     44    double lattice_distortion, double radius, 
    4545    double sld, double solvent_sld) 
    4646{ 
    … …  
    6767            const double qa = qab*cos_phi; 
    6868            const double qb = qab*sin_phi; 
    69             const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     69            const double form = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
    7070            inner_sum += GAUSS_W[j] * form; 
    7171        } 
    … …  
    7777    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    7878 
    79     return sc_volume_fraction(radius, dnn) * Pq * Zq; 
     79    return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    8080} 
    8181 
    8282static double 
    8383Iqabc(double qa, double qb, double qc, 
    84     double dnn, double d_factor, double radius, 
     84    double lattice_spacing, double lattice_distortion, double radius, 
    8585    double sld, double solvent_sld) 
    8686{ 
    8787    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
    8888    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    89     const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
    90     return sc_volume_fraction(radius, dnn) * Pq * Zq; 
     89    const double Zq = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 
     90    return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 
    9191} 

  • sasmodels/models/sc_paracrystal.py

    rda7b26b rb3f4831  
    11r""" 
    2 .. warning:: This model and this model description are under review following  
    3              concerns raised by SasView users. If you need to use this model,  
    4              please email help@sasview.org for the latest situation. *The  
     2.. warning:: This model and this model description are under review following 
     3             concerns raised by SasView users. If you need to use this model, 
     4             please email help@sasview.org for the latest situation. *The 
    55             SasView Developers. September 2018.* 
    6               
     6 
    77Definition 
    88---------- 
    … …  
    104104 
    105105Authorship and Verification 
    106 --------------------------- 
     106---------------------------- 
    107107 
    108108* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    109 * **Last Modified by:** Paul Butler **Date:** September 29, 2016 
    110 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 
     109* **Last Modified by:** Paul Butler **Date:** September 16, 2018 
     110* **Last Reviewed by:** Paul Butler **Date:** September 16, 2018 
    111111""" 
    112112 
    … …  
    129129        scale: volume fraction of spheres 
    130130        bkg:background, R: radius of sphere 
    131         dnn: Nearest neighbor distance 
    132         d_factor: Paracrystal distortion factor 
     131        lattice_spacing: Nearest neighbor distance 
     132        lattice_distortion: Paracrystal distortion factor 
    133133        radius: radius of the spheres 
    134134        sldSph: SLD of the sphere 
    … …  
    139139# pylint: disable=bad-whitespace, line-too-long 
    140140#             ["name", "units", default, [lower, upper], "type","description"], 
    141 parameters = [["dnn",         "Ang",       220.0, [0.0, inf],  "",            "Nearest neighbor distance"], 
    142               ["d_factor",    "",           0.06, [-inf, inf], "",            "Paracrystal distortion factor"], 
     141parameters = [["lattice_spacing",         "Ang",       220.0, [0.0, inf],  "",  "Lattice spacing"], 
     142              ["lattice_distortion",    "",           0.06, [-inf, inf], "",   "Paracrystal distortion factor"], 
    143143              ["radius",      "Ang",        40.0, [0.0, inf],  "volume",      "Radius of sphere"], 
    144144              ["sld",  "1e-6/Ang^2",         3.0, [0.0, inf],  "sld",         "Sphere scattering length density"], 
    … …  
    155155    # copied from bcc_paracrystal 
    156156    radius = 10**np.random.uniform(1.3, 4) 
    157     d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
    158     dnn_fraction = np.random.beta(a=10, b=1) 
    159     dnn = radius*4/np.sqrt(4)/dnn_fraction 
     157    lattice_distortion = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
     158    lattice_spacing_fraction = np.random.beta(a=10, b=1) 
     159    lattice_spacing = radius*4/np.sqrt(4)/lattice_spacing_fraction 
    160160    pars = dict( 
    161161        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
    162         dnn=dnn, 
    163         d_factor=d_factor, 
     162        lattice_spacing=lattice_spacing, 
     163        lattice_distortion=lattice_distortion, 
    164164        radius=radius, 
    165165    ) 
Note: See TracChangeset for help on using the changeset viewer.