Changeset 3a45c2c in sasmodels


Ignore:
Timestamp:
Nov 20, 2017 11:50:21 AM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
costrafo411
Parents:
65ceb7d (diff), fa70e04 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into costrafo411

Files:
11 added
2 deleted
26 edited

Legend:

Unmodified
Added
Removed

  • doc/genmodel.py

    r745b7bb rb866abf  
    33import sys, os, math, re 
    44import numpy as np 
     5import matplotlib 
     6matplotlib.use('Agg') 
    57import matplotlib.pyplot as plt 
    68sys.path.insert(0, os.path.abspath('..')) 

  • doc/guide/index.rst

    r2e66ef5 rc0d7ab3  
    99   intro.rst 
    1010   install.rst 
     11   gpu_setup.rst 
    1112   pd/polydispersity.rst 
    1213   resolution.rst 

  • doc/guide/install.rst

    rf8a2baa rc0d7ab3  
    5353This will allow you to edit the files in the package and have the changes 
    5454show up immediately in python the next time you load your program. 
    55  
    56 OpenCL Installation 
    57 ******************* 
    58 *Warning! GPU devices do not in general offer the same level of memory 
    59 protection as CPU devices. If your code attempts to write outside allocated 
    60 memory buffers unpredicatable behaviour may result (eg, your video display 
    61 may freeze, or your system may crash, etc). Do not install OpenCL drivers 
    62 without first checking for known issues (eg, some computer manufacturers 
    63 install modified graphics drivers so replacing these may not be a good 
    64 idea!). If in doubt, seek advice from an IT professional before proceeding 
    65 further.* 
    66  
    67 Check if you have OpenCL already installed 
    68 ========================================== 
    69  
    70 **Windows** 
    71  
    72 The following instructions are based on 
    73 http://web.engr.oregonstate.edu/~mjb/cs475/DoIHaveOpenCL.pdf 
    74  
    75 * Go to: Start -> Control Panel -> System & Security -> Administrative Tools 
    76 * Double Click on Computer Managment 
    77 * Click on Device Manager 
    78 * Click open Display Adapters 
    79 * Right-click on available adapter and select Properties 
    80 * Click on Driver 
    81 * Go to Driver Details 
    82 * Scroll down and see if OpenCL is installed (look for OpenCL*.dll files) 
    83  
    84 **Mac OSX** 
    85  
    86 For OS X operating systems higher than 10.6 OpenCL is shipped along with 
    87 the system. 
    88  
    89 However, OpenCL has had a rocky history on Macs. Apple provide a useful 
    90 compatibility table at https://support.apple.com/en-us/HT202823 
    91  
    92  
    93 Installation 
    94 ============ 
    95  
    96 **Windows** 
    97  
    98 Depending on the graphic card in your system, drivers 
    99 can be obtained from different sources: 
    100  
    101 * NVIDIA: https://developer.nvidia.com/opencl 
    102 * AMD: http://developer.amd.com/tools-and-sdks/opencl-zone/ 
    103  
    104  
    105 **Mac OSX** 
    106  
    107 N/A 
    108  
    109 You cannot download OpenCL driver updates for your Mac. They are packaged 
    110 with the normal quarterly OS X updates from Apple. 
    111  
    112  
    113 .. note:: 
    114     Intel provides OpenCL drivers for Intel processors at 
    115     https://software.intel.com/en-us/articles/opencl-drivers 
    116     These can sometimes make use of special vector instructions across multiple 
    117     processors, so it is worth installing if the GPU does not support double 
    118     precision. You can install this driver alongside the GPU driver for NVIDIA 
    119     or AMD. 
    120  
    121  
    122 GPU Selection 
    123 ************* 
    124  
    125 sasmodels evaluations can run on your graphics card (GPU) or they can run 
    126 on the processor (CPU). In general, calculations performed on the GPU will run faster. 
    127  
    128 To run on the GPU, your computer must have OpenCL drivers installed. 
    129 For information about OpenCL installation see this 
    130 :ref:`opencl-installation` guidance. 
    131  
    132 Where the model is evaluated is a little bit complicated. 
    133 If the model has the line *single=False* then it requires double precision. 
    134 If the GPU is single precision only, then it will try running via OpenCL 
    135 on the CPU.  If the OpenCL driver is not available for the CPU then 
    136 it will run as a normal program on the CPU. 
    137  
    138 For models with a large number of parameters or with a lot of code, 
    139 the GPU may be too small to run the program effectively. 
    140 In this case, you should try simplifying the model, maybe breaking it 
    141 into several different modules so that you don't need *IF* statements in your code. 
    142 If it is still too big, you can set *opencl=False* in the model file and 
    143 the model will only run as a normal program on the CPU. 
    144 This will not usually be necessary. 
    145  
    146 Device Selection 
    147 ================ 
    148 If you have multiple GPU devices you can tell SasView which device to use. 
    149 By default, SasView looks for one GPU and one CPU device 
    150 from available OpenCL platforms. 
    151  
    152 SasView prefers AMD or NVIDIA drivers for GPU, and prefers Intel or 
    153 Apple drivers for CPU. Both GPU and CPU are included on the assumption that CPU 
    154 is always available and supports double precision. 
    155  
    156 The device order is important: GPU is checked before CPU on the assumption that 
    157 it will be faster. By examining ~/sasview.log you can see which device SasView 
    158 chose to run the model. 
    159  
    160 **If you don't want to use OpenCL, you can set** *SAS_OPENCL=None* 
    161 **in your environment settings, and it will only use normal programs.** 
    162  
    163 If you want to use one of the other devices, you can run the following 
    164 from the python console in SasView:: 
    165  
    166     import pyopencl as cl 
    167     cl.create_some_context() 
    168  
    169 This will provide a menu of different OpenCL drivers available. 
    170 When one is selected, it will say "set PYOPENCL_CTX=..." 
    171 Use that value as the value of *SAS_OPENCL*. 
    172  
    173 Compiler Selection 
    174 ================== 
    175 For models run as normal programs, you may need to specify a compiler. 
    176 This is done using the SAS_COMPILER environment variable. 
    177 Set it to *tinycc* for the tinycc compiler, *msvc* for the 
    178 Microsoft Visual C compiler, or *mingw* for the MinGW compiler. 
    179 TinyCC is provided with SasView so that is the default. 
    180 If you want one of the other compilers, be sure to have it available 
    181 in your *PATH* so SasView can find it! 
    182  
    183  
    184 *Document History* 
    185  
    186 | 2017-05-17 Paul Kienzle 

  • doc/guide/magnetism/magnetism.rst

    r64eecf7 r1f058ea  
    7777===========   ================================================================ 
    7878 M0_sld        = $D_M M_0$ 
    79  Up_theta      = $\theta_{up}$ 
     79 Up_theta      = $\theta_\mathrm{up}$ 
    8080 M_theta       = $\theta_M$ 
    8181 M_phi         = $\phi_M$ 
    … …  
    8787    The values of the 'Up_frac_i' and 'Up_frac_f' must be in the range 0 to 1. 
    8888 
    89 .. note:: 
    90     This help document was last changed by Steve King, 02May2015 
     89*Document History* 
    9190 
    92 * Document History * 
    93  
     91| 2015-05-02 Steve King 
    9492| 2017-05-08 Paul Kienzle 

  • doc/guide/pd/polydispersity.rst

    rf8a2baa r1f058ea  
    9595           \exp\left(-\frac{(x - \bar x)^2}{2\sigma^2}\right) 
    9696 
    97 where $\bar x$ is the mean of the distribution and *Norm* is a normalization factor 
    98 which is determined during the numerical calculation. 
     97where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
     98factor which is determined during the numerical calculation. 
    9999 
    100100The polydispersity is 
    … …  
    122122during the numerical calculation. 
    123123 
    124 The median value for the distribution will be the value given for the respective 
    125 size parameter, for example, *radius=60*. 
     124The median value for the distribution will be the value given for the 
     125respective size parameter, for example, *radius=60*. 
    126126 
    127127The polydispersity is given by $\sigma$ 
    … …  
    208208 
    209209Many commercial Dynamic Light Scattering (DLS) instruments produce a size 
    210 polydispersity parameter, sometimes even given the symbol $p$ This 
     210polydispersity parameter, sometimes even given the symbol $p$\ ! This 
    211211parameter is defined as the relative standard deviation coefficient of 
    212212variation of the size distribution and is NOT the same as the polydispersity 

  • doc/guide/plugin.rst

    r30b60d2 r3048ec6  
    7878     at the start of the model documentation and as a tooltip in the SasView GUI 
    7979 
    80 - a **short discription**: 
     80- a **short description**: 
    8181   - this is the **description** string in the *.py* file 
    8282   - this is a medium length description which appears when you click 
    … …  
    233233 
    234234**name = "mymodel"** defines the name of the model that is shown to the user. 
    235 If it is not provided, it will use the name of the model file, with '_' 
    236 replaced by spaces and the parts capitalized.  So *adsorbed_layer.py* will 
    237 become *Adsorbed Layer*.  The predefined models all use the name of the 
    238 model file as the name of the model, so the default may be changed. 
     235If it is not provided it will use the name of the model file. The name must 
     236be a valid variable name, starting with a letter and contains only letters, 
     237numbers or underscore.  Spaces, dashes, and other symbols are not permitted. 
    239238 
    240239**title = "short description"** is short description of the model which 
    … …  
    298297- **"name"** is the name of the parameter shown on the FitPage. 
    299298 
     299  - the name must be a valid variable name, starting with a letter and 
     300    containing only letters, numbers and underscore. 
     301 
    300302  - parameter names should follow the mathematical convention; e.g., 
    301303    *radius_core* not *core_radius*, or *sld_solvent* not *solvent_sld*. 
    … …  
    366368    dispersion. 
    367369 
     370Some models will have integer parameters, such as number of pearls in the 
     371pearl necklace model, or number of shells in the multi-layer vesicle model. 
     372The optimizers in BUMPS treat all parameters as floating point numbers which 
     373can take arbitrary values, even for integer parameters, so your model should 
     374round the incoming parameter value to the nearest integer inside your model 
     375you should round to the nearest integer.  In C code, you can do this using:: 
     376 
     377    static double 
     378    Iq(double q, ..., double fp_n, ...) 
     379    { 
     380        int n = (int)(fp_n + 0.5); 
     381        ... 
     382    } 
     383 
     384in python:: 
     385 
     386    def Iq(q, ..., n, ...): 
     387        n = int(n+0.5) 
     388        ... 
     389 
     390Derivative based optimizers such as Levenberg-Marquardt will not work 
     391for integer parameters since the partial derivative is always zero, but 
     392the remaining optimizers (DREAM, differential evolution, Nelder-Mead simplex) 
     393will still function. 
    368394 
    369395Model Computation 
    … …  
    391417     \approx \frac{\sum_{i=1}^n G(x_i) I(q,x_i)}{\sum_{i=1}^n G(x_i) V(x_i)} 
    392418 
    393 That is, the indivdual models do not need to include polydispersity 
     419That is, the individual models do not need to include polydispersity 
    394420calculations, but instead rely on numerical integration to compute the 
    395421appropriately smeared pattern.   Angular dispersion values over polar angle 
    … …  
    411437The parameters *par1, par2, ...* are the list of non-orientation parameters 
    412438to the model in the order that they appear in the parameter table. 
    413 **Note that the autogenerated model file uses** *x* **rather than** *q*. 
     439**Note that the auto-generated model file uses** *x* **rather than** *q*. 
    414440 
    415441The *.py* file should import trigonometric and exponential functions from 
    … …  
    935961across multiple parameters can be very slow). 
    936962 
    937 If your model has 2D orientational calculation, then you should also 
     963If your model has 2D orientation calculation, then you should also 
    938964test with:: 
    939965 

  • doc/guide/resolution.rst

    r30b60d2 r1f058ea  
    1717resolution contribution into a model calculation/simulation (which by definition 
    1818will be exact) to make it more representative of what has been measured 
    19 experimentally - a process called *smearing*. sasmodels does the latter. 
     19experimentally - a process called *smearing*. Sasmodels does the latter. 
    2020 
    2121Both smearing and desmearing rely on functions to describe the resolution 
    22 effect. sasmodels provides three smearing algorithms: 
     22effect. Sasmodels provides three smearing algorithms: 
    2323 
    2424*  *Slit Smearing* 
    … …  
    9999 
    100100For discrete $q$ values, at the $q$ values of the data points and at the $q$ 
    101 values extended up to $q_N = q_i + \Delta q_v$ the smeared 
     101values extended up to $q_N = q_i + \Delta q_u$ the smeared 
    102102intensity can be approximately calculated as 
    103103 

  • doc/guide/scripting.rst

    r2e66ef5 r4aa5dce  
    6969    $ bumps example/model.py --preview 
    7070 
     71Note that bumps and sasmodels are included as part of the SasView 
     72distribution.  On windows, bumps can be called from the cmd prompt 
     73as follows:: 
     74 
     75    SasViewCom bumps.cli example/model.py --preview 
     76 
    7177Using sasmodels directly 
    7278======================== 
    … …  
    105111    plt.loglog(q, Iq) 
    106112    plt.show() 
     113 
     114On windows, this can be called from the cmd prompt using sasview as:: 
     115 
     116    SasViewCom example/cylinder_eval.py 

  • example/oriented_usans.py

    r1cd24b4 r74b0495  
    66 
    77# Spherical particle data, not ellipsoids 
    8 sans, usans = load_data('../../sasview/sasview/test/1d_data/latex_smeared.xml') 
     8sans, usans = load_data('latex_smeared.xml', index='all') 
    99usans.qmin, usans.qmax = np.min(usans.x), np.max(usans.x) 
    1010usans.mask = (usans.x < 0.0) 

  • example/simul_fit.py

    r1a4d4c0 r74b0495  
    1 #!/usr/bin/env python 
    2 # -*- coding: utf-8 -*- 
    3  
    4 # To Sasview/documents/scripts 
    5  
    61from bumps.names import * 
    72from sasmodels.core import load_model 
    … …  
    94from sasmodels.data import load_data, plot_data 
    105 
     6# latex data, same sample usans and sans 
     7# particles radius ~2300, uniform dispersity 
     8datasets = load_data('latex_smeared.xml', index='all') 
     9#[print(data) for data in datasets] 
    1110 
    12 """ IMPORT THE DATA USED """ 
    13 datafiles = ['latex_smeared_out_0.txt', 'latex_smeared_out_1.txt'] 
    14 datasets = [load_data(el) for el in datafiles] 
     11# A single sphere model to share between the datasets.  We will use 
     12# FreeVariables below to set the parameters that are independent between 
     13# the datasets. 
     14kernel = load_model('sphere') 
     15pars = dict(scale=0.01, background=0.0, sld=5.0, sld_solvent=0.0, radius=1500., 
     16            #radius_pd=0.1, radius_pd_n=35, 
     17            ) 
     18model = Model(kernel, **pars) 
    1519 
    16 for data in datasets: 
    17     data.qmin = 0.0 
    18     data.qmax = 10.0 
     20# radius and polydispersity (if any) are shared 
     21model.radius.range(0, inf) 
     22#model.radius_pd.range(0, 1) 
    1923 
    20 #sphere model 
    21 kernel = load_model('sphere', dtype="single") 
    22 pars = dict(scale=0.01, background=0.0, sld=1.0, sld_solvent=6.0, radius=1500.) 
    23 model = Model(kernel, **pars) 
    24 model.radius.range(0, inf) 
    25 #model.background.range(-inf, inf) 
    26 #model.scale.range(0, inf) 
    27 model.sld.range(-inf, inf) 
    28 model.sld_solvent.range(-inf, inf) 
     24# Contrast and dilution are the same for both measurements, but are not 
     25# separable with a single measurement (i.e., I(q) ~ F(q) contrast^2 Vf), 
     26# so fit one of scale, sld or solvent sld.  With absolute scaling from 
     27# data reduction, can use the same parameter for both datasets. 
     28model.scale.range(0, inf) 
     29#model.sld.range(-inf, inf) 
     30#model.sld_solvent.range(-inf, inf) 
    2931 
     32# Background is different for sans and usans so set it as a free variable 
     33# in the model. 
    3034free = FreeVariables( 
    31     names=[data.filename for data in datasets], 
     35    names=[data.run[0] for data in datasets], 
    3236    background=model.background, 
    33     scale=model.scale, 
    3437    ) 
    3538free.background.range(-inf, inf) 
    36 free.scale.range(0, inf) 
    3739 
    38 M = [Experiment(data=data, model=model) for data in datasets] 
     40# Note: can access the parameters for the individual models using 
     41# free.background[0] and free.background[1], setting constraints or 
     42# ranges as appropriate. 
     43 
     44# For more complex systems where different datasets require independent models, 
     45# separate models can be defined, with parameters tied together using 
     46# constraint expressions.  For example, the following could be used to fit 
     47# data set 1 to spheres and data set 2 to cylinders of the same volume: 
     48#    model1 = Model(load_model('sphere')) 
     49#    model2 = Model(load_model('cylinder')) 
     50#    model1.sld = model2.sld 
     51#    model1.sld_solvent = model2.sld_solvent 
     52#    model1.scale = model2.scale 
     53#    # set cylinders and spheres to the same volume 
     54#    model1.radius = (3/4*model2.radius**2*model2.length)**(1/3) 
     55#    model1.background.range(0, 2) 
     56#    model2.background.range(0, 2) 
     57 
     58# Setup the experiments, sharing the same model across all datasets. 
     59M = [Experiment(data=data, model=model, name=data.run[0]) for data in datasets] 
    3960 
    4061problem = FitProblem(M, freevars=free) 
    41  
    42 print(problem._parameters) 

  • explore/precision.py

    r487e695 r237c9cf  
    8181        elif xrange == "linear": 
    8282            lin_min, lin_max, lin_steps = 1, 1000, 2000 
     83            lin_min, lin_max, lin_steps = 0.001, 2, 2000 
    8384        elif xrange == "log": 
    8485            log_min, log_max, log_steps = -3, 5, 400 
    … …  
    321322    np_function=lambda x: np.fmod(x, 2*np.pi), 
    322323    ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"), 
     324) 
     325add_function( 
     326    name="debye", 
     327    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
     328    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
     329    ocl_function=make_ocl(""" 
     330    const double qsq = q*q; 
     331    if (qsq < 1.0) { // Pade approximation 
     332        const double x = qsq; 
     333        if (0) { // 0.36 single 
     334            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}] 
     335            return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1); 
     336        } else if (0) { // 1.0 for single 
     337            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}] 
     338            const double A1=1./24., A2=1./84, A3=-1./3360; 
     339            const double B1=3./8., B2=3./56., B3=1./336.; 
     340            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
     341        } else if (1) { // 1.0 for single, 0.25 for double 
     342            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     343            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     344            const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; 
     345            return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) 
     346                  /((((B4*x + B3)*x + B2)*x + B1)*x + 1.); 
     347        } else { // 1.0 for single, 0.5 for double 
     348            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     349            const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.; 
     350            const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.; 
     351            return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.) 
     352                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
     353        } 
     354    } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     355        const double x = qsq; 
     356        const double C0 = +1.; 
     357        const double C1 = -1./3.; 
     358        const double C2 = +1./12.; 
     359        const double C3 = -1./60.; 
     360        const double C4 = +1./360.; 
     361        const double C5 = -1./2520.; 
     362        const double C6 = +1./20160.; 
     363        const double C7 = -1./181440.; 
     364        //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
     365        //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
     366        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
     367    } else { 
     368        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); 
     369    } 
     370    """, "sas_debye"), 
    323371) 
    324372 

  • sasmodels/bumps_model.py

    r3330bb4 r74b0495  
    133133    """ 
    134134    _cache = None # type: Dict[str, np.ndarray] 
    135     def __init__(self, data, model, cutoff=1e-5): 
     135    def __init__(self, data, model, cutoff=1e-5, name=None): 
    136136        # type: (Data, Model, float) -> None 
    137137        # remember inputs so we can inspect from outside 
     138        self.name = data.filename if name is None else name 
    138139        self.model = model 
    139140        self.cutoff = cutoff 
    … …  
    204205        """ 
    205206        data, theory, resid = self._data, self.theory(), self.residuals() 
    206         plot_theory(data, theory, resid, view, Iq_calc=self.Iq_calc) 
     207        # TODO: hack to display oriented usans 2-D pattern 
     208        Iq_calc = self.Iq_calc if isinstance(self.Iq_calc, tuple) else None 
     209        plot_theory(data, theory, resid, view, Iq_calc=Iq_calc) 
    207210 
    208211    def simulate_data(self, noise=None): 

  • sasmodels/compare.py

    r765eb0e rced5bd2  
    7373    -res=0 sets the resolution width dQ/Q if calculating with resolution 
    7474    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0 
     75    -q=min:max alternative specification of qrange 
    7576    -nq=128 sets the number of Q points in the data set 
    7677    -1d*/-2d computes 1d or 2d data 
    … …  
    533534                % (pd, n, nsigma, nsigma, pdtype) 
    534535    if M0 != 0.: 
    535         line += "  M0:%.3f  mphi:%.1f  mtheta:%.1f" % (M0, mphi, mtheta) 
     536        line += "  M0:%.3f  mtheta:%.1f  mphi:%.1f" % (M0, mtheta, mphi) 
    536537    return line 
    537538 
    … …  
    768769    'accuracy', 'is2d' and 'view' parsed from the command line. 
    769770    """ 
    770     qmax, nq, res = opts['qmax'], opts['nq'], opts['res'] 
     771    qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 
    771772    if opts['is2d']: 
    772773        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray 
    … …  
    777778    else: 
    778779        if opts['view'] == 'log' and not opts['zero']: 
    779             qmax = math.log10(qmax) 
    780             q = np.logspace(qmax-3, qmax, nq) 
     780            q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 
    781781        else: 
    782             q = np.linspace(0.001*qmax, qmax, nq) 
     782            q = np.linspace(qmin, qmax, nq) 
    783783        if opts['zero']: 
    784784            q = np.hstack((0, q)) 
    … …  
    955955        else: 
    956956            err, errstr, errview = abs(relerr), "rel err", "log" 
     957            if (err == 0.).all(): 
     958                errview = 'linear' 
    957959        if 0:  # 95% cutoff 
    958960            sorted = np.sort(err.flatten()) 
    … …  
    960962            err[err > cutoff] = cutoff 
    961963        #err,errstr = base/comp,"ratio" 
    962         plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    963         if view == 'linear': 
    964             plt.xscale('linear') 
     964        plot_theory(data, None, resid=err, view=view, use_data=use_data) 
     965        plt.yscale(errview) 
    965966        plt.title("max %s = %.3g"%(errstr, abs(err).max())) 
    966967        #cbar_title = errstr if errview=="linear" else "log "+errstr 
    … …  
    10011002 
    10021003    # Data generation 
    1003     'data=', 'noise=', 'res=', 
    1004     'nq=', 'lowq', 'midq', 'highq', 'exq', 'zero', 
     1004    'data=', 'noise=', 'res=', 'nq=', 'q=', 
     1005    'lowq', 'midq', 'highq', 'exq', 'zero', 
    10051006    '2d', '1d', 
    10061007 
    … …  
    11221123        'view'      : 'log', 
    11231124        'is2d'      : False, 
     1125        'qmin'      : None, 
    11241126        'qmax'      : 0.05, 
    11251127        'nq'        : 128, 
    … …  
    11591161        elif arg == '-zero':    opts['zero'] = True 
    11601162        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:]) 
     1163        elif arg.startswith('-q='): 
     1164            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 
    11611165        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
    11621166        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    … …  
    12071211 
    12081212    # Create the computational engines 
     1213    if opts['qmin'] is None: 
     1214        opts['qmin'] = 0.001*opts['qmax'] 
    12091215    if opts['datafile'] is not None: 
    12101216        data = load_data(os.path.expanduser(opts['datafile'])) 

  • sasmodels/conversion_table.py

    r0c2da4b r505d0ad  
    1111account for that. 
    1212 
    13 Usage: 
    14 <old_Sasview_version> : { 
    15     <new_model_name> : [ 
    16         <old_model_name> , 
    17         { 
    18             <new_param_name_1> : <old_param_name_1>, 
    19             ... 
    20             <new_param_name_n> : <old_param_name_n> 
    21         } 
    22     ] 
    23 } 
     13Usage:: 
     14 
     15    <old_Sasview_version> : { 
     16        <new_model_name> : [ 
     17            <old_model_name> , 
     18            { 
     19                <new_param_name_1> : <old_param_name_1>, 
     20                ... 
     21                <new_param_name_n> : <old_param_name_n> 
     22            } 
     23        ] 
     24    } 
    2425 
    2526Any future parameter and model name changes can and should be given in this 

  • sasmodels/custom/__init__.py

    r3c852d4 r0f48f1e  
    3131        if fullname in sys.modules: 
    3232            del sys.modules[fullname] 
     33        if path.endswith(".py") and os.path.exists(path) and os.path.exists(path+"c"): 
     34            # remove automatic pyc file before loading a py file 
     35            os.unlink(path+"c") 
    3336        module = imp.load_source(fullname, os.path.expanduser(path)) 
    34         #os.unlink(path+"c")  # remove the automatic pyc file 
    3537        return module 
    3638 

  • sasmodels/data.py

    r630156b rba7302a  
    4444    Data = Union["Data1D", "Data2D", "SesansData"] 
    4545 
    46 def load_data(filename): 
     46def load_data(filename, index=0): 
    4747    # type: (str) -> Data 
    4848    """ 
    … …  
    5555        filename, indexstr = filename[:-1].split('[') 
    5656        index = int(indexstr) 
    57     else: 
    58         index = None 
    5957    datasets = loader.load(filename) 
    60     if datasets is None: 
     58    if not datasets:  # None or [] 
    6159        raise IOError("Data %r could not be loaded" % filename) 
    6260    if not isinstance(datasets, list): 
    6361        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 
    73     return data 
     62    for data in datasets: 
     63        if hasattr(data, 'x'): 
     64            data.qmin, data.qmax = data.x.min(), data.x.max() 
     65            data.mask = (np.isnan(data.y) if data.y is not None 
     66                        else np.zeros_like(data.x, dtype='bool')) 
     67        elif hasattr(data, 'qx_data'): 
     68            data.mask = ~data.mask 
     69    return datasets[index] if index != 'all' else datasets 
    7470 
    7571 
    … …  
    438434 
    439435    scale = data.x**4 if view == 'q4' else 1.0 
     436    xscale = yscale = 'linear' if view == 'linear' else 'log' 
    440437 
    441438    if use_data or use_theory: 
    … …  
    470467            plt.ylim(*limits) 
    471468 
    472         plt.xscale('linear' if not some_present or non_positive_x 
    473                    else view if view is not None 
    474                    else 'log') 
    475         plt.yscale('linear' 
    476                    if view == 'q4' or not some_present or not all_positive 
    477                    else view if view is not None 
    478                    else 'log') 
     469 
     470        xscale = ('linear' if not some_present or non_positive_x 
     471                  else view if view is not None 
     472                  else 'log') 
     473        yscale = ('linear' 
     474                  if view == 'q4' or not some_present or not all_positive 
     475                  else view if view is not None 
     476                  else 'log') 
     477        plt.xscale(xscale) 
    479478        plt.xlabel("$q$/A$^{-1}$") 
     479        plt.yscale(yscale) 
    480480        plt.ylabel('$I(q)$') 
    481481        title = ("data and model" if use_theory and use_data 
    … …  
    505505        plt.xlabel("$q$/A$^{-1}$") 
    506506        plt.ylabel('residuals') 
    507         plt.xscale('linear') 
    508507        plt.title('(model - Iq)/dIq') 
     508        plt.xscale(xscale) 
     509        plt.yscale('linear') 
    509510 
    510511 

  • sasmodels/generate.py

    r30b60d2 r6db17bd  
    370370    """ 
    371371    # Note: need 0xffffffff&val to force an unsigned 32-bit number 
     372    try: 
     373        source = source.encode('utf8') 
     374    except AttributeError: # bytes has no encode attribute in python 3 
     375        pass 
    372376    return "%08X"%(0xffffffff&crc32(source)) 
    373377 

  • sasmodels/kernelpy.py

    r883ecf4 r3b32f8d  
    3434        self.info = model_info 
    3535        self.dtype = np.dtype('d') 
     36        logging.info("load python model " + self.info.name) 
    3637 
    3738    def make_kernel(self, q_vectors): 
    38         logging.info("creating python kernel " + self.info.name) 
    3939        q_input = PyInput(q_vectors, dtype=F64) 
    4040        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c rc69d6d6  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    … …  
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //return length_a * length_b * length_c; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
    … …  
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     37 
    3738    const double mu = 0.5 * q * length_b; 
    38      
     39 
    3940    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     41    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     42    //double vol = length_a * length_b * length_c + 
     43    //       2.0 * thick_rim_a * length_b * length_c + 
    4344    //       2.0 * thick_rim_b * length_a * length_c + 
    4445    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     46 
    4647    // Scale sides by B 
    4748    const double a_scaled = length_a / length_b; 
    4849    const double c_scaled = length_c / length_b; 
    4950 
    50     // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 
    51     // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 
    52     // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 
    53     // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 
    54     double ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    55     double tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
     51    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
     52    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
     53    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
    5654 
    5755    double Vin = length_a * length_b * length_c; 
    … …  
    6260    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    6361    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
     62    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
     63    double Vot = Vin + V1 + V2 + V3; 
    6464 
    6565    // Scale factors (note that drC is not used later) 
    … …  
    6767    const double drhoA = (arim_sld-solvent_sld); 
    6868    const double drhoB = (brim_sld-solvent_sld); 
    69     //const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     69    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     70 
    7071 
    7172    // Precompute scale factors for combining cross terms from the shape 
    7273    const double scale23 = drhoA*V1; 
    7374    const double scale14 = drhoB*V2; 
    74     const double scale12 = drho0*Vin - scale23 - scale14; 
     75    const double scale24 = drhoC*V3; 
     76    const double scale11 = drho0*Vin; 
     77    const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
    7578 
    7679    // outer integral (with gauss points), integration limits = 0, 1 
    … …  
    8386        // inner integral (with gauss points), integration limits = 0, 1 
    8487        double inner_total = 0.0; 
     88        double inner_total_crim = 0.0; 
    8589        for(int j=0; j<76; j++) { 
    8690            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    … …  
    8892            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    8993            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    90             const double si2 = sas_sinx_x(mu_proj * cos_uu); 
     94            const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    9195            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    9296            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
    … …  
    9498            // Expression in libCylinder.c (neither drC nor Vot are used) 
    9599            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
     100            const double form_crim = scale11*si1*si2; 
    96101 
    97             // To note also that in csparallelepiped.cpp, there is a function called 
    98             // cspkernel, which seems to make more sense and has the following comment: 
    99             //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010. 
    100             //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)* 
    101             //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    102             // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     102 
    105103            //  correct FF : sum of square of phase factors 
    106104            inner_total += Gauss76Wt[j] * form * form; 
     105            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
    107106        } 
    108107        inner_total *= 0.5; 
    109  
     108        inner_total_crim *= 0.5; 
    110109        // now sum up the outer integral 
    111110        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    112         outer_total += Gauss76Wt[i] * inner_total * si * si; 
     111        const double si_crim = sas_sinx_x(mu * tc * sigma); 
     112        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
    113113    } 
    114114    outer_total *= 0.5; 
    … …  
    154154 
    155155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    156     // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 
    157     // but for the moment I let it like this until understanding better the code. 
     156    // the scaling by B. 
    158157    double ta = length_a + 2.0*thick_rim_a; 
    159     double tb = length_a + 2.0*thick_rim_b; 
    160     double tc = length_a + 2.0*thick_rim_c; 
     158    double tb = length_b + 2.0*thick_rim_b; 
     159    double tc = length_c + 2.0*thick_rim_c; 
    161160    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162161    double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    … …  
    166165    double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167166    double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     167 
    169168 
    170169    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    173172               + drA*(siAt-siA)*siB*siC*V1 
    174173               + drB*siA*(siBt-siB)*siC*V2 
    175                + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     174               + drC*siA*siB*(siCt-siC)*V3); 
     175 
    177176    return 1.0e-4 * f * f; 
    178177} 

  • sasmodels/models/core_shell_parallelepiped.py

    r8f04da4 rfa70e04  
    211211 
    212212# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
    213 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    214 tests = [[{}, 0.2, 0.533149288477], 
    215          [{}, [0.2], [0.533149288477]], 
    216          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
    217          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
    218         ] 
    219 del qx, qy  # not necessary to delete, but cleaner 
     213if 0:  # pak: model rewrite; need to update tests 
     214    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     215    tests = [[{}, 0.2, 0.533149288477], 
     216            [{}, [0.2], [0.533149288477]], 
     217            [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
     218            [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
     219            ] 
     220    del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/flexible_cylinder.py

    r31df0c9 r2573fa1  
    9090    length = 10**np.random.uniform(2, 6) 
    9191    radius = 10**np.random.uniform(1, 3) 
    92     kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments 
     92    kuhn_length = 10**np.random.uniform(-2, 0)*length 
    9393    pars = dict( 
    9494        length=length, 
    … …  
    100100tests = [ 
    101101    # Accuracy tests based on content in test/utest_other_models.py 
    102     # Currently fails in OCL 
    103     # [{'length':     1000.0, 
    104     #  'kuhn_length': 100.0, 
    105     #  'radius':       20.0, 
    106     #  'sld':           1.0, 
    107     #  'sld_solvent':   6.3, 
    108     #  'background':    0.0001, 
    109     #  }, 0.001, 3509.2187], 
     102    [{'length':     1000.0,  # test T1 
     103      'kuhn_length': 100.0, 
     104      'radius':       20.0, 
     105      'sld':           1.0, 
     106      'sld_solvent':   6.3, 
     107      'background':    0.0001, 
     108     }, 0.001, 3509.2187], 
    110109 
    111110    # Additional tests with larger range of parameters 
    112     [{'length':    1000.0, 
     111    [{'length':    1000.0,  # test T2 
    113112      'kuhn_length': 100.0, 
    114113      'radius':       20.0, 
    … …  
    117116      'background':    0.0001, 
    118117     }, 1.0, 0.000595345], 
    119     [{'length':        10.0, 
     118    [{'length':        10.0,  # test T3 
    120119      'kuhn_length': 800.0, 
    121120      'radius':        2.0, 
    … …  
    124123      'background':    0.001, 
    125124     }, 0.1, 1.55228], 
    126     [{'length':        100.0, 
     125    [{'length':        100.0,  # test T4 
    127126      'kuhn_length': 800.0, 
    128127      'radius':       50.0, 
    … …  
    133132    ] 
    134133 
     134# There are a few branches in the code that ought to have test values: 
     135# 
     136# For length > 4 * kuhn_length 
     137#        if length > 10 * kuhn_length then C is scaled by 3.06 (L/b)^(-0.44) 
     138#        q*kuhn_length <= 3.1  => Sexv_new 
     139#           dS/dQ < 0 has different behaviour from dS/dQ >= 0 
     140#  T2    q*kuhn_length > 3.1   => a_long 
     141# 
     142# For length <= 4 * kuhn_length 
     143#        q*kuhn_length <= max(1.9/Rg_short, 3.0)  => Sdebye((q*Rg)^2) 
     144#           q*Rg < 0.5 uses Pade approx, q*Rg > 1.0 uses math lib 
     145#  T3,T4 q*kuhn_length > max(1.9/Rg_short, 3.0)   => a_short 
     146# 
     147# Note that the transitions between branches may be abrupt.  You can see a 
     148# several percent change around length=10*kuhn_length and length=4*kuhn_length 
     149# using the following: 
     150# 
     151#    sascomp flexible_cylinder -calc=double -sets=10 length=10*kuhn_length,10.000001*kuhn_length 
     152#    sascomp flexible_cylinder -calc=double -sets=10 length=4*kuhn_length,4.000001*kuhn_length 
     153# 
     154# The transition between low q and high q around q*kuhn_length = 3 seems 
     155# to be good to 4 digits or better.  This was tested by computing the value 
     156# on each branches near the transition point and reporting the relative error 
     157# for kuhn lengths of 10, 100 and 1000 and a variety of length:kuhn_length 
     158# ratios. 

  • sasmodels/models/lib/wrc_cyl.c

    r92ce163 r18a2bfc  
    22    Functions for WRC implementation of flexible cylinders 
    33*/ 
    4 double Sk_WR(double q, double L, double b); 
    5  
    6  
    7 static double 
    8 AlphaSquare(double x) 
    9 { 
    10     // Potentially faster. Needs proper benchmarking. 
    11     // add native_powr to kernel_template 
    12     //double t = native_powr( (1.0 + (x/3.12)*(x/3.12) + 
    13     //     (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
    14     //return t; 
    15  
    16     return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0); 
    17 } 
    18  
    19 // 
    20 static double 
    21 Rgsquarezero(double q, double L, double b) 
    22 { 
    23     const double r = b/L; 
    24     return (L*b/6.0) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 
    25  
    26 } 
    27  
    28 // 
    29 static double 
    30 Rgsquareshort(double q, double L, double b) 
    31 { 
    32     return AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
    33 } 
    34  
    35 // 
    36 static double 
    37 Rgsquare(double q, double L, double b) 
    38 { 
    39     return AlphaSquare(L/b)*L*b/6.0; 
    40 } 
    41  
    42 static double 
    43 sech_WR(double x) 
    44 { 
    45     return(1/cosh(x)); 
    46 } 
    47  
    48 static double 
    49 a1long(double q, double L, double b, double p1, double p2, double q0) 
    50 { 
    51     double C; 
    52  
    53     if( L/b > 10.0) { 
    54         C = 3.06/pow((L/b),0.44); 
     4 
     5static double 
     6Rgsquare(double L, double b) 
     7{ 
     8    const double x = L/b; 
     9    // Use Horner's method to evaluate Pedersen eq 15: 
     10    //     alpha^2 = [1.0 + (x/3.12)^2 + (x/8.67)^3] ^ (0.176/3) 
     11    const double alphasq = 
     12        pow(1.0 + x*x*(1.027284681130835e-01 + 1.534414548417740e-03*x), 
     13            5.866666666666667e-02); 
     14    return alphasq*L*b/6.0; 
     15} 
     16 
     17static double 
     18Rgsquareshort(double L, double b) 
     19{ 
     20    const double r = b/L;  // = 1/n_b in Pedersen ref. 
     21    return Rgsquare(L, b) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 
     22} 
     23 
     24static double 
     25w_WR(double x) 
     26{ 
     27    // Pedersen eq. 16: 
     28    //    w = [1 + tanh((x-C4)/C5)]/2 
     29    const double C4 = 1.523; 
     30    const double C5 = 0.1477; 
     31    return 0.5 + 0.5*tanh((x - C4)/C5); 
     32} 
     33 
     34static double 
     35Sdebye(double qsq) 
     36{ 
     37#if FLOAT_SIZE>4 
     38#  define DEBYE_CUTOFF 0.25  // 1e-15 error 
     39#else 
     40#  define DEBYE_CUTOFF 1.0  // 4e-7 error 
     41#endif 
     42 
     43    if (qsq < DEBYE_CUTOFF) { 
     44        const double x = qsq; 
     45        // mathematica: PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     46        const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     47        const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; 
     48        return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) 
     49             / ((((B4*x + B3)*x + B2)*x + B1)*x + 1.); 
    5550    } else { 
    56         C = 1.0; 
     51        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); 
    5752    } 
    58  
     53} 
     54 
     55static double 
     56a_long(double q, double L, double b/*, double p1, double p2, double q0*/) 
     57{ 
     58    const double p1 = 4.12; 
     59    const double p2 = 4.42; 
     60    const double q0 = 3.1; 
     61 
     62    // Constants C1, ..., C5 derived from least squares fit (Pedersen, eq 13,16) 
    5963    const double C1 = 1.22; 
    6064    const double C2 = 0.4288; 
    … …  
    6468    const double miu = 0.585; 
    6569 
    66     const double Rg2 = Rgsquare(q,L,b); 
    67     const double Rg22 = Rg2*Rg2; 
    68     const double Rg = sqrt(Rg2); 
    69     const double Rgb = Rg*q0/b; 
    70  
    71     const double b2 = b*b; 
     70    const double C = (L/b>10.0 ? 3.06*pow(L/b, -0.44) : 1.0); 
     71    //printf("branch B-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 
     72    const double r2 = Rgsquare(L,b); 
     73    const double r = sqrt(r2); 
     74    const double qr_b = q0*r/b; 
     75    const double qr_b_sq = qr_b*qr_b; 
     76    const double qr_b_4 = qr_b_sq*qr_b_sq; 
     77    const double qr_b_miu = pow(qr_b, -1.0/miu); 
     78    const double em1_qr_b_sq = expm1(-qr_b_sq); 
     79    const double sech2 = 1.0/square(cosh((qr_b-C4)/C5)); 
     80    const double w = w_WR(qr_b); 
     81 
     82    const double t1 = pow(q0, 1.0 + p1 + p2)/(b*(p1-p2)); 
     83    const double t2 = C/(15.0*L) * ( 
     84        + 14.0*b*b*em1_qr_b_sq/(q0*qr_b_sq) 
     85        + 2.0*q0*r2*exp(-qr_b_sq)*(11.0 + 7.0/qr_b_sq)); 
     86    const double t11 = ((C3*qr_b_miu + C2)*qr_b_miu + C1)*qr_b_miu; 
     87    const double t3 = r*sech2/(2.*C5)*t11; 
     88    const double t4 = r*(em1_qr_b_sq + qr_b_sq)*sech2 / (C5*qr_b_4); 
     89    const double t5 = -4.0*r*qr_b*em1_qr_b_sq/qr_b_4 * (1.0 - w); 
     90    const double t10 = 2.0*(em1_qr_b_sq + qr_b_sq)/qr_b_4 * (1.0 - w); //=Sdebye*(1-w) 
     91    const double t6 = 4.0*b/q0 * t10; 
     92    const double t7 = r*((-3.0*C3*qr_b_miu -2.0*C2)*qr_b_miu -1.0*C1)*qr_b_miu/(miu*qr_b); 
     93    const double t9 = C*b/L * (4.0 - exp(-qr_b_sq) * (11.0 + 7.0/qr_b_sq) + 7.0/qr_b_sq)/15.0; 
     94    const double t12 = b*b*M_PI/(L*q0*q0) + t2 + t3 - t4 + t5 - t6 + t7*w; 
     95    const double t13 = -b*M_PI/(L*q0) + t9 + t10 + t11*w; 
     96 
     97    const double a1 = pow(q0,p1)*t13 - t1*pow(q0,-p2)*(t12 + b*p1/q0*t13); 
     98    const double a2 = t1*pow(q0,-p1)*(t12 + b*p1/q0*t13); 
     99 
     100    const double ans = a1*pow(q*b, -p1) + a2*pow(q*b, -p2) + M_PI/(q*L); 
     101    return ans; 
     102} 
     103 
     104static double 
     105_short(double r2, double exp_qr_b, double L, double b, double p1short, 
     106        double p2short, double q0) 
     107{ 
     108    const double qr2 = q0*q0 * r2; 
    72109    const double b3 = b*b*b; 
    73     const double b4 = b3*b; 
    74     const double q02 = q0*q0; 
    75     const double q03 = q0*q0*q0; 
    76     const double q04 = q03*q0; 
    77     const double q05 = q04*q0; 
    78  
    79     const double Rg02 = Rg2*q02; 
    80  
    81     const double t1 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2))) * 
    82          ((11.0/15.0 + (7.0*b2)/(15.0*Rg02))) + 
    83          (7.0*b2)/(15.0*Rg02)))); 
    84  
    85     const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    86          Rg02/b2))*((1.0 + 0.5*(((-1.0) - 
    87          tanh((-C4 + Rgb/C5))))))); 
    88  
    89     const double t3 = ((C3*pow(Rgb,((-3.0)/miu)) + 
    90          C2*pow(Rgb,((-2.0)/miu)) + 
    91          C1*pow(Rgb,((-1.0)/miu)))); 
    92  
    93     const double t4 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
    94  
    95     const double t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - 
    96          b*p2*pow(q0,((-1.0) - p1 - p2)))); 
    97  
    98     const double t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + 
    99          (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 
    100          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 
    101          (7.0*b2)/(15.0*Rg02)))*Rg2)/b))); 
    102  
    103     const double t7 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    104          C2*pow(((Rgb)),((-2.0)/miu)) + 
    105          C1*pow(((Rgb)),((-1.0)/miu))))*pow(sech_WR(((-C4) + 
    106          Rgb)/C5),2)); 
    107  
    108     const double t8 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    109          Rg02/b2))*pow(sech_WR(((-C4) + Rgb)/C5),2)); 
    110  
    111     const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    112          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) - 
    113          tanh(((-C4) + Rgb)/C5)))))); 
    114  
    115     const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    116          Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    117          Rgb)/C5)))))); 
    118  
    119     const double t11 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 
    120           3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),((-1.0) - 
    121           2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),((-1.0) - 
    122           1.0/miu)))/miu)); 
    123  
    124     const double t12 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
    125  
    126     const double t13 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 
    127           (7.0*b2)/(15.0*q02* Rg2))) + 
    128           (7.0*b2)/(15.0*Rg02)))); 
    129  
    130     const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    131           Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    132           Rgb)/C5)))))); 
    133  
    134     const double t15 = ((C3*pow(((Rgb)),((-3.0)/miu)) + 
    135         C2*pow(((Rgb)),((-2.0)/miu)) + 
    136         C1*pow(((Rgb)),((-1.0)/miu)))); 
    137  
    138  
    139     double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 
    140         0.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
    141         (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 
    142         t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) - 
    143         b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 
    144         t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) + 
    145         Rgb)/C5))))))))))); 
    146  
    147     return (yy); 
    148 } 
    149  
    150 static double 
    151 a2long(double q, double L, double b, double p1, double p2, double q0) 
    152 { 
    153     double C; 
    154  
    155     if( L/b > 10.0) { 
    156         C = 3.06/pow((L/b),0.44); 
    157     } else { 
    158         C = 1.0; 
    159     } 
    160  
    161     const double C1 = 1.22; 
    162     const double C2 = 0.4288; 
    163     const double C3 = -1.651; 
    164     const double C4 = 1.523; 
    165     const double C5 = 0.1477; 
    166     const double miu = 0.585; 
    167  
    168     const double Rg2 = Rgsquare(q,L,b); 
    169     const double Rg22 = Rg2*Rg2; 
    170     const double b2 = b*b; 
    171     const double b3 = b*b*b; 
    172     const double b4 = b3*b; 
    173     const double q02 = q0*q0; 
    174     const double q03 = q0*q0*q0; 
    175     const double q04 = q03*q0; 
    176     const double q05 = q04*q0; 
    177     const double Rg = sqrt(Rg2); 
    178     const double Rgb = Rg*q0/b; 
    179     const double Rg02 = Rg2*q02; 
    180  
    181     const double t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - 
    182          b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
    183  
    184     const double t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + 
    185          (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 
    186          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 
    187          (7*b2)/(15.0*Rg02)))*Rg2)/b)))/L; 
    188  
    189     const double t3 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    190          C2*pow(((Rgb)),((-2.0)/miu)) + 
    191          C1*pow(((Rgb)),((-1.0)/miu))))* 
    192          pow(sech_WR(((-C4) +Rgb)/C5),2.0))/(2.0*C5); 
    193  
    194     const double t4 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    195          Rg02/b2))*pow(sech_WR(((-C4) + 
    196          Rgb)/C5),2))/(C5*q04*Rg22); 
    197  
    198     const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    199          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))* 
    200          ((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    201          Rgb)/C5))))))/(q04*Rg22); 
    202  
    203     const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    204          Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    205          Rgb)/C5))))))/(q05*Rg22); 
    206  
    207     const double t7 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 
    208          3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)), 
    209          ((-1.0) - 2.0/miu)))/miu - (C1*Rg*pow(((Rgb)), 
    210          ((-1.0) - 1.0/miu)))/miu)); 
    211  
    212     const double t8 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
    213  
    214     const double t9 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 
    215          (7.0*b2)/(15*Rg02))) + (7.0*b2)/(15.0*Rg02))))/L; 
    216  
    217     const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) + 
    218           Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    219           Rgb)/C5))))))/(q04*Rg22); 
    220  
    221     const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) + 
    222          t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
    223          (((-((b*M_PI)/(L*q0))) + t9 + t10 + 
    224          0.5*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    225          C2*pow(((Rgb)),((-2.0)/miu)) + 
    226          C1*pow(((Rgb)),((-1.0)/miu))))* 
    227          ((1.0 + tanh(((-C4) + Rgb)/C5)))))))))); 
    228  
    229     return (yy); 
    230 } 
    231  
    232 // 
    233 static double 
    234 a_short(double q, double L, double b, double p1short, 
    235         double p2short, double factor, double pdiff, double q0) 
    236 { 
    237     const double Rg2_sh = Rgsquareshort(q,L,b); 
    238     const double Rg2_sh2 = Rg2_sh*Rg2_sh; 
    239     const double b3 = b*b*b; 
    240     const double t1 = ((q0*q0*Rg2_sh)/(b*b)); 
    241     const double Et1 = exp(t1); 
    242     const double Emt1 = 1.0/Et1; 
    243     const double q02 = q0*q0; 
    244     const double q0p = pow(q0,(-4.0 + p1short) ); 
    245  
    246     double yy = ((factor/(L*(pdiff)*Rg2_sh2)* 
    247         ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 
    248         2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 
    249         2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*M_PI*q02*q0*Rg2_sh2 + 
    250         Et1*p2short*M_PI*q02*q0*Rg2_sh2)))))); 
    251  
    252     return(yy); 
    253 } 
    254 static double 
    255 a1short(double q, double L, double b, double p1short, double p2short, double q0) 
    256 { 
    257  
    258     double factor = 1.0; 
    259     return a_short(q, L, b, p1short, p2short, factor, p1short - p2short, q0); 
    260 } 
    261  
    262 static double 
    263 a2short(double q, double L, double b, double p1short, double p2short, double q0) 
    264 { 
    265     double factor = -1.0; 
    266     return a_short(q, L, b, p2short, p1short, factor, p1short-p2short, q0); 
    267 } 
    268  
    269 //WR named this w (too generic) 
    270 static double 
    271 w_WR(double x) 
    272 { 
    273     return 0.5*(1 + tanh((x - 1.523)/0.1477)); 
    274 } 
    275  
    276 // 
    277 static double 
    278 u1(double q, double L, double b) 
    279 { 
    280     return Rgsquareshort(q,L,b)*q*q; 
    281 } 
    282  
    283 static double 
    284 u_WR(double q, double L, double b) 
    285 { 
    286     return Rgsquare(q,L,b)*q*q; 
    287 } 
    288  
    289 static double 
    290 Sdebye_kernel(double arg) 
    291 { 
    292     // ORIGINAL 
    293     double result = 2.0*(exp(-arg) + arg -1.0)/(arg*arg); 
    294  
    295     // CONVERSION 1 from http://herbie.uwplse.org/ 
    296     // 
    297     // exhibits discontinuity - needs more investigation 
    298     //double a1 = 1.0/6.0; 
    299     //double a2 = 1.0/72.0; 
    300     //double a3 = 1.0/24.0; 
    301     //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2); 
    302  
    303     return result; 
    304 } 
    305 static double 
    306 Sdebye(double q, double L, double b) 
    307 { 
    308     double arg = u_WR(q,L,b); 
    309     return Sdebye_kernel(arg); 
    310 } 
    311  
    312 // 
    313 static double 
    314 Sdebye1(double q, double L, double b) 
    315 { 
    316     double arg = u1(q,L,b); 
    317     return Sdebye_kernel(arg); 
    318  
    319 } 
    320  
    321 // 
     110    const double q0p = pow(q0, -4.0 + p1short); 
     111 
     112    double yy = 1.0/(L*r2*r2) * b/exp_qr_b*q0p 
     113        * (8.0*b3*L 
     114           - 8.0*b3*exp_qr_b*L 
     115           + 2.0*b3*exp_qr_b*L*p2short 
     116           - 2.0*b*exp_qr_b*L*p2short*qr2 
     117           + 4.0*b*exp_qr_b*L*qr2 
     118           - 2.0*b3*L*p2short 
     119           + 4.0*b*L*qr2 
     120           - M_PI*exp_qr_b*qr2*q0*r2 
     121           + M_PI*exp_qr_b*p2short*qr2*q0*r2); 
     122 
     123    return yy; 
     124} 
     125static double 
     126a_short(double qp, double L, double b 
     127        /*double p1short, double p2short*/, double q0) 
     128{ 
     129    const double p1short = 5.36; 
     130    const double p2short = 5.62; 
     131 
     132    const double r2 = Rgsquareshort(L,b); 
     133    const double exp_qr_b = exp(r2*square(q0/b)); 
     134    const double pdiff = p1short - p2short; 
     135    const double a1 = _short(r2,exp_qr_b,L,b,p1short,p2short,q0)/pdiff; 
     136    const double a2= -_short(r2,exp_qr_b,L,b,p2short,p1short,q0)/pdiff; 
     137    const double ans = a1*pow(qp*b, -p1short) + a2*pow(qp*b, -p2short) + M_PI/(qp*L); 
     138    return ans; 
     139} 
     140 
     141 
    322142static double 
    323143Sexv(double q, double L, double b) 
    324144{ 
    325  
     145    // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w 
    326146    const double C1=1.22; 
    327147    const double C2=0.4288; 
    328148    const double C3=-1.651; 
    329149    const double miu = 0.585; 
    330     const double qRg = q*sqrt(Rgsquare(q,L,b)); 
    331     const double x = pow(qRg, -1.0/miu); 
    332  
    333     double yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) + 
    334             w_WR(qRg)*x*(C1 + x*(C2 + x*C3)); 
    335     return (yy); 
    336 } 
    337  
    338  
    339 static double 
    340 Sexvnew(double q, double L, double b) 
    341 { 
    342     double yy; 
    343  
    344     const double C1 =1.22; 
    345     const double C2 =0.4288; 
    346     const double C3 =-1.651; 
    347     const double miu = 0.585; 
     150    const double qr = q*sqrt(Rgsquare(L,b)); 
     151    const double qr_miu = pow(qr, -1.0/miu); 
     152    const double w = w_WR(qr); 
     153    const double t10 = Sdebye(qr*qr)*(1.0 - w); 
     154    const double t11 = ((C3*qr_miu + C2)*qr_miu + C1)*qr_miu; 
     155 
     156    return t10 + w*t11; 
     157} 
     158 
     159// Modified by Yun on Oct. 15, 
     160static double 
     161Sexv_new(double q, double L, double b) 
     162{ 
     163    const double qr = q*sqrt(Rgsquare(L,b)); 
     164    const double qr2 = qr*qr; 
     165    const double C = (L/b > 10.0) ? 3.06*pow(L/b, -0.44) : 1.0; 
     166    const double t9 = C*b/L * (4.0 - exp(-qr2) * (11.0 + 7.0/qr2) + 7.0/qr2)/15.0; 
     167 
     168    const double Sexv_orig = Sexv(q, L, b); 
     169 
     170    // calculating the derivative to decide on the correction (cutoff) term? 
     171    // Note: this is modified from WRs original code 
    348172    const double del=1.05; 
    349     const double qRg = q*sqrt(Rgsquare(q,L,b)); 
    350     const double x = pow(qRg, -1.0/miu); 
    351  
    352  
    353     //calculating the derivative to decide on the corection (cutoff) term? 
    354     // I have modified this from WRs original code 
    355     const double qdel = (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q); 
    356     const double C_star2 =(qdel >= 0.0) ? 0.0 : 1.0; 
    357  
    358     yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) + 
    359          C_star2*w_WR(qRg)* 
    360          x*(C1 + x*(C2 + x*C3)); 
    361  
    362     return (yy); 
    363 } 
    364  
    365 double Sk_WR(double q, double L, double b) 
    366 { 
    367     // 
    368     const double p1 = 4.12; 
    369     const double p2 = 4.42; 
    370     const double p1short = 5.36; 
    371     const double p2short = 5.62; 
    372     const double q0 = 3.1; 
    373  
    374     double q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
    375     double Sexvmodify, ans; 
    376  
    377     const double C =(L/b > 10.0) ? 3.06/pow((L/b),0.44) : 1.0; 
    378  
    379     if( L > 4*b ) { // Longer Chains 
    380        if (q*b <= 3.1) {   //Modified by Yun on Oct. 15, 
    381          Sexvmodify = Sexvnew(q, L, b); 
    382          ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - 
    383             (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
    384  
    385        } else { //q(i)*b > 3.1 
    386          ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + 
    387                a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + M_PI/(q*L); 
    388        } 
    389     } else { //L <= 4*b Shorter Chains 
    390        if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
    391          if (q*b<=0.01) { 
    392             ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
    393          } else { 
    394             ans = Sdebye1(q,L,b); 
    395          } 
    396        } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
    397          ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + 
    398                a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + 
    399                M_PI/(q*L); 
    400        } 
     173    const double qdel = (Sexv(q*del,L,b) - Sexv_orig)/(q*(del - 1.0)); 
     174 
     175    if (qdel < 0) { 
     176        //printf("branch A1-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 
     177        return t9 + Sexv_orig; 
     178    } else { 
     179        //printf("branch A2-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 
     180        const double w = w_WR(qr); 
     181        const double t10 = Sdebye(qr*qr)*(1.0 - w); 
     182        return t9 + t10; 
    401183    } 
    402  
    403   return(ans); 
    404 } 
     184} 
     185 
     186 
     187static double 
     188Sk_WR(double q, double L, double b) 
     189{ 
     190    const double Rg_short = sqrt(Rgsquareshort(L, b)); 
     191    double q0short = fmax(1.9/Rg_short, 3.0); 
     192    double ans; 
     193 
     194 
     195    if( L > 4*b ) { // L > 4*b : Longer Chains 
     196        if (q*b <= 3.1) { 
     197            ans = Sexv_new(q, L, b); 
     198        } else { //q(i)*b > 3.1 
     199            ans = a_long(q, L, b /*, p1, p2, q0*/); 
     200        } 
     201    } else { // L <= 4*b : Shorter Chains 
     202        if (q*b <= q0short) { // q*b <= fmax(1.9/Rg_short, 3) 
     203            //printf("branch C-%d q=%g L=%g b=%g\n", square(q*Rg_short)<DEBYE_CUTOFF, q, L, b); 
     204            // Note that q0short is usually 3, but it will be greater than 3 
     205            // small enough b, depending on the L/b ratio: 
     206            //     L/b == 1 => b < 2.37 
     207            //     L/b == 2 => b < 1.36 
     208            //     L/b == 3 => b < 1.00 
     209            //     L/b == 4 => b < 0.816 
     210            // 2017-10-01 pkienzle: moved low q approximation into Sdebye() 
     211            ans = Sdebye(square(q*Rg_short)); 
     212        } else {  // q*b > max(1.9/Rg_short, 3) 
     213            //printf("branch D q=%g L=%g b=%g\n", q, L, b); 
     214            ans = a_short(q, L, b /*, p1short, p2short*/, q0short); 
     215        } 
     216    } 
     217 
     218    return ans; 
     219} 

  • sasmodels/product.py

    r058460c r146793b  
    100100    # Remember the component info blocks so we can build the model 
    101101    model_info.composition = ('product', [p_info, s_info]) 
     102    model_info.control = p_info.control 
     103    model_info.hidden = p_info.hidden 
     104    if getattr(p_info, 'profile', None) is not None: 
     105        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters) 
     106        def profile(**kwargs): 
     107            # extract the profile args 
     108            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars) 
     109            return p_info.profile(**kwargs) 
     110    else: 
     111        profile = None 
     112    model_info.profile = profile 
     113    model_info.profile_axes = p_info.profile_axes 
     114 
    102115    # TODO: delegate random to p_info, s_info 
    103116    #model_info.random = lambda: {} 
    … …  
    129142    def __init__(self, model_info, P, S): 
    130143        # type: (ModelInfo, KernelModel, KernelModel) -> None 
     144        #: Combined info plock for the product model 
    131145        self.info = model_info 
     146        #: Form factor modelling individual particles. 
    132147        self.P = P 
     148        #: Structure factor modelling interaction between particles. 
    133149        self.S = S 
    134         self.dtype = P.dtype 
     150        #: Model precision. This is not really relevant, since it is the 
     151        #: individual P and S models that control the effective dtype, 
     152        #: converting the q-vectors to the correct type when the kernels 
     153        #: for each are created. Ideally this should be set to the more 
     154        #: precise type to avoid loss of precision, but precision in q is 
     155        #: not critical (single is good enough for our purposes), so it just 
     156        #: uses the precision of the form factor. 
     157        self.dtype = P.dtype  # type: np.dtype 
    135158 
    136159    def make_kernel(self, q_vectors): 

  • sasmodels/sasview_model.py

    r9644b5a redb0f85  
    1616import logging 
    1717from os.path import basename, splitext, abspath, getmtime 
    18 import thread 
     18try: 
     19    import _thread as thread 
     20except ImportError: 
     21    import thread 
    1922 
    2023import numpy as np  # type: ignore 
    … …  
    202205                                           structure_factor._model_info) 
    203206    ConstructedModel = make_model_from_info(model_info) 
    204     return ConstructedModel() 
     207    return ConstructedModel(form_factor.multiplicity) 
    205208 
    206209 
    … …  
    320323    #: True if model has multiplicity 
    321324    is_multiplicity_model = False 
    322     #: Mulitplicity information 
     325    #: Multiplicity information 
    323326    multiplicity_info = None # type: MultiplicityInfoType 
    324327 
    … …  
    351354        # and lines to plot. 
    352355 
    353         # Get the list of hidden parameters given the mulitplicity 
     356        # Get the list of hidden parameters given the multiplicity 
    354357        # Don't include multiplicity in the list of parameters 
    355358        self.multiplicity = multiplicity 

  • sasmodels/direct_model.py

    rd1ff3a5 r65ceb7d  
    202202 
    203203        if self.data_type == 'sesans': 
    204             q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
    205             index = slice(None, None) 
    206             res = None 
    207             if data.y is not None: 
    208                 Iq, dIq = data.y, data.dy 
     204            from sas.sascalc.data_util.nxsunit import Converter 
     205            qmax, qunits = data.sample.zacceptance 
     206            SElength = Converter(data._xunit)(data.x, "A") 
     207            zaccept = Converter(qunits)(qmax, "1/A"), 
     208            Rmax = 10000000 
     209            index = slice(None, None)  # equivalent to index [:] 
     210            Iq = data.y[index] 
     211            dIq = data.dy[index] 
     212            oriented = getattr(data, 'oriented', False) 
     213            if oriented: 
     214                res = sesans.OrientedSesansTransform(data.x[index], SElength, zaccept, Rmax) 
     215                # Oriented sesans transform produces q_calc = [qx, qy] 
     216                q_vectors = res.q_calc 
    209217            else: 
    210                 Iq, dIq = None, None 
    211             #self._theory = np.zeros_like(q) 
    212             q_vectors = [q] 
    213             q_mono = sesans.make_all_q(data) 
     218                res = sesans.SesansTransform(data.x[index], SElength, zaccept, Rmax) 
     219                # Unoriented sesans transform produces q_calc = q 
     220                q_vectors = [res.q_calc] 
    214221        elif self.data_type == 'Iqxy': 
    215222            #if not model.info.parameters.has_2d: 
    … …  
    230237            #self._theory = np.zeros_like(self.Iq) 
    231238            q_vectors = res.q_calc 
    232             q_mono = [] 
    233239        elif self.data_type == 'Iq': 
    234240            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    … …  
    255261            #self._theory = np.zeros_like(self.Iq) 
    256262            q_vectors = [res.q_calc] 
    257             q_mono = [] 
    258263        elif self.data_type == 'Iq-oriented': 
    259264            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    … …  
    272277                                      qy_width=data.dxl[index]) 
    273278            q_vectors = res.q_calc 
    274             q_mono = [] 
    275279        else: 
    276280            raise ValueError("Unknown data type") # never gets here 
    … …  
    279283        # so we can save/restore state 
    280284        self._kernel_inputs = q_vectors 
    281         self._kernel_mono_inputs = q_mono 
    282285        self._kernel = None 
    283286        self.Iq, self.dIq, self.index = Iq, dIq, index 
    … …  
    317320        if self._kernel is None: 
    318321            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    319             self._kernel_mono = ( 
    320                 self._model.make_kernel(self._kernel_mono_inputs) 
    321                 if self._kernel_mono_inputs else None) 
    322322 
    323323        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    … …  
    326326        # TODO: extend plotting of calculate Iq to other measurement types 
    327327        # TODO: refactor so we don't store the result in the model 
    328         self.Iq_calc = Iq_calc 
    329         if self.data_type == 'sesans': 
    330             Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 
    331                        if self._kernel_mono_inputs else None) 
    332             result = sesans.transform(self._data, 
    333                                       self._kernel_inputs[0], Iq_calc, 
    334                                       self._kernel_mono_inputs, Iq_mono) 
    335         else: 
    336             result = self.resolution.apply(Iq_calc) 
    337             if hasattr(self.resolution, 'nx'): 
    338                 self.Iq_calc = ( 
    339                     self.resolution.qx_calc, self.resolution.qy_calc, 
    340                     np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
    341                 ) 
     328        self.Iq_calc = None 
     329        result = self.resolution.apply(Iq_calc) 
     330        if hasattr(self.resolution, 'nx'): 
     331            self.Iq_calc = ( 
     332                self.resolution.qx_calc, self.resolution.qy_calc, 
     333                np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
     334            ) 
    342335        return result 
    343336 

  • sasmodels/sesans.py

    r9f91afe r65ceb7d  
    2020    Spin-Echo SANS transform calculator.  Similar to a resolution function, 
    2121    the SesansTransform object takes I(q) for the set of *q_calc* values and 
    22     produces a transformed dataset 
     22    produces a transformed dataset. 
    2323 
    2424    *SElength* (A) is the set of spin-echo lengths in the measured data. 
    … …  
    4848 
    4949    def apply(self, Iq): 
    50         # tye: (np.ndarray) -> np.ndarray 
    5150        G0 = np.dot(self._H0, Iq) 
    5251        G = np.dot(self._H.T, Iq) 
    … …  
    7877        self.q_calc = q 
    7978        self._H, self._H0 = H, H0 
     79 
     80class OrientedSesansTransform(object): 
     81    """ 
     82    Oriented Spin-Echo SANS transform calculator.  Similar to a resolution 
     83    function, the OrientedSesansTransform object takes I(q) for the set 
     84    of *q_calc* values and produces a transformed dataset. 
     85 
     86    *SElength* (A) is the set of spin-echo lengths in the measured data. 
     87 
     88    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 
     89    echo encoding dimension (for ToF: Q of min(R) and max(lam)). 
     90 
     91    *Rmax* (A) is the maximum size sensitivity; larger radius requires more 
     92    computation time. 
     93    """ 
     94    #: SElength from the data in the original data units; not used by transform 
     95    #: but the GUI uses it, so make sure that it is present. 
     96    q = None  # type: np.ndarray 
     97 
     98    #: q values to calculate when computing transform 
     99    q_calc = None  # type: np.ndarray 
     100 
     101    # transform arrays 
     102    _cosmat = None  # type: np.ndarray 
     103    _cos0 = None # type: np.ndarray 
     104    _Iq_shape = None # type: Tuple[int, int] 
     105 
     106    def __init__(self, z, SElength, zaccept, Rmax): 
     107        # type: (np.ndarray, float, float) -> None 
     108        #import logging; logging.info("creating SESANS transform") 
     109        self.q = z 
     110        self._set_cosmat(SElength, zaccept, Rmax) 
     111 
     112    def apply(self, Iq): 
     113        dq = self.q_calc[0][0] 
     114        Iq = np.reshape(Iq, self._Iq_shape) 
     115        G0 = self._cos0 * np.sum(Iq) * dq 
     116        G = np.sum(np.dot(Iq, self._cosmat), axis=0) * dq 
     117        P = G - G0 
     118        return P 
     119 
     120    def _set_cosmat(self, SElength, zaccept, Rmax): 
     121        # type: (np.ndarray, float, float) -> None 
     122        # Force float32 arrays, otherwise run into memory problems on some machines 
     123        SElength = np.asarray(SElength, dtype='float32') 
     124 
     125        # Rmax = #value in text box somewhere in FitPage? 
     126        q_max = 2 * pi / (SElength[1] - SElength[0]) 
     127        q_min = 0.1 * 2 * pi / (np.size(SElength) * SElength[-1]) 
     128        q_min *= 100 
     129 
     130        q = np.arange(q_min, q_max, q_min, dtype='float32') 
     131        dq = q_min 
     132 
     133        cos0 = np.float32(dq / (2 * pi)) 
     134        cosmat = np.float32(dq / (2 * pi)) * np.cos(q[:, None] * SElength[None, :]) 
     135 
     136        qx, qy = np.meshgrid(q, q) 
     137        self._Iq_shape = qx.shape 
     138        self.q_calc = qx.flatten(), qy.flatten() 
     139        self._cosmat, self._cos0 = cosmat, cos0 
Note: See TracChangeset for help on using the changeset viewer.