Changeset 7904790 in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 8:40:15 AM (8 years ago)
Author:
wojciech
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3b12dea
Parents:
a5af4e1 (diff), 52a3db3 (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' of ​https://github.com/SasView/sasmodels

Files:
4 added
20 edited
7 moved

Legend:

Unmodified
Added
Removed

  • doc/genmodel.py

    rc094758 r044a9c1  
    22import numpy as np 
    33import matplotlib.pyplot as plt 
     4import pylab 
    45sys.path.insert(0, os.path.abspath('..')) 
    56from sasmodels import generate, core 
    … …  
    2829    'q_max'     : 1.0, 
    2930    'nq'        : 1000, 
    30     'nq2d'      : 400, 
     31    'nq2d'      : 100, 
    3132    'vmin'      : 1e-3,  # floor for the 2D data results 
    3233    'qx_max'    : 0.5, 
    … …  
    5960    if opts['zscale'] == 'log': 
    6061        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 
    61     h = ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='upper', 
    62            extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=ice_cm()) 
    63     # , vmin=vmin, vmax=vmax) 
     62    ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='lower', 
     63        extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=pylab.cm.jet)    
    6464    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 
    6565    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 
     66 
     67#            im = self.subplot.imshow(output, interpolation='nearest', 
     68#                                     origin='lower', 
     69#                                     vmin=zmin_temp, vmax=self.zmax_2D, 
     70#                                     cmap=self.cmap, 
     71#                                     extent=(self.xmin_2D, self.xmax_2D, 
     72#                                             self.ymin_2D, self.ymax_2D)) 
    6673 
    6774def ice_cm(): 
    … …  
    7784 
    7885 
    79 # Generate image (comment IF for 1D/2D for the moment) and generate only 1D 
     86# Generate image  
    8087fig_height = 3.0 # in 
    8188fig_left = 0.6 # in 
    … …  
    97104    plot_2d(model, opts, ax2d) 
    98105    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     106    plot_1d(model, opts, ax1d) 
    99107    #ax.set_aspect('square') 
    100108else: 
    … …  
    109117    fig = plt.figure(figsize=aspect) 
    110118    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
    111 plot_1d(model, opts, ax1d) 
     119    plot_1d(model, opts, ax1d) 
    112120 
    113121# Save image in model/img 
    … …  
    121129captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n' 
    122130captionstr += '\n' 
    123 captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
     131if model_info['has_2d']: 
     132    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n' 
     133else: 
     134    captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    124135captionstr += '\n' 
    125136 

  • example/sesansfit.py

    r84db7a5 r170ea69  
    11from bumps.names import * 
    22from sasmodels import core, bumps_model, sesans 
     3from sas.sascalc.dataloader.loader import Loader 
    34 
    45HAS_CONVERTER = True 
    … …  
    89    HAS_CONVERTER = False 
    910 
     11 
    1012def get_bumps_model(model_name): 
    1113    kernel = core.load_model(model_name) 
    … …  
    1315    return model 
    1416 
    15 def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[]): 
     17def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[], acceptance_angle=None): 
    1618    """ 
     19 
    1720    @param file: SESANS file location 
    1821    @param model: Bumps model object or model name - can be model, model_1 * model_2, and/or model_1 + model_2 
    … …  
    2427    """ 
    2528    try: 
    26         from sas.sascalc.dataloader.loader import Loader 
    2729        loader = Loader() 
    2830        data = loader.load(file) 
    … …  
    5557            dy = err_data 
    5658            sample = Sample() 
     59            needs_all_q = acceptance_angle is not None 
    5760        data = SESANSData1D() 
     61        data.acceptance_angle = acceptance_angle 
    5862 
     63    data.needs_all_q = acceptance_angle is not None 
    5964    if "radius" in initial_vals: 
    6065        radius = initial_vals.get("radius") 

  • example/sesansfit.sh

    rbdb653c r170ea69  
    11#!/bin/bash 
    22 
    3 SASVIEW=$PWD/../../sasview/src 
    4 PYTHONPATH=$PWD/../:$PWD/../../bumps:$PWD/../../periodictable:$SASVIEW 
     3# Need to fix the paths to sasmodels and sasview if no eggs present 
     4echo $PWD 
     5ONEUP="$(dirname "$PWD")" 
     6PROJECTS="$(dirname "$ONEUP")" 
     7CCOLON="C:/" 
     8CSLASH="/c/" 
     9SASMODELSBASE=$PROJECTS/sasmodels/ 
     10SASMODELS="${SASMODELSBASE/$CSLASH/$CCOLON}" 
     11SASVIEWBASE=$PROJECTS/sasview/src/ 
     12SASVIEW="${SASVIEWBASE/$CSLASH/$CCOLON}" 
     13PYTHONPATH="$SASVIEW;$SASMODELS" 
    514export PYOPENCL_CTX PYTHONPATH 
    615 

  • sasmodels/convert.py

    r78d3341 r74f7238  
    1515    'be_polyelectrolyte', 
    1616    'correlation_length', 
    17     'fractal_core_shell' 
     17    'fractal_core_shell', 
    1818    'binary_hard_sphere', 
    1919] 

  • sasmodels/core.py

    r7b3e62c rd5ba841  
    170170    """ 
    171171    value, weight = zip(*pars) 
    172     if len(value) > 1: 
    173         value = [v.flatten() for v in np.meshgrid(*value)] 
    174         weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)]) 
    175         weight = np.prod(weight, axis=0) 
     172    value = [v.flatten() for v in np.meshgrid(*value)] 
     173    weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)]) 
     174    weight = np.prod(weight, axis=0) 
    176175    return value, weight 
    177176 
    178 def call_kernel(kernel, pars, cutoff=0): 
     177def call_kernel(kernel, pars, cutoff=0, mono=False): 
    179178    """ 
    180179    Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
    … …  
    189188    fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    190189                  for name in kernel.fixed_pars] 
    191     pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
     190    if mono: 
     191        pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
     192                   for name in kernel.pd_pars] 
     193    else: 
     194        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
    192195    return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
    193196 

  • sasmodels/direct_model.py

    r17bbadd r02e70ff  
    6969 
    7070        if self.data_type == 'sesans': 
     71             
    7172            q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
    7273            index = slice(None, None) 
    … …  
    7778                Iq, dIq = None, None 
    7879            #self._theory = np.zeros_like(q) 
    79             q_vectors = [q] 
     80            q_vectors = [q]             
     81            q_mono = sesans.make_all_q(data) 
    8082        elif self.data_type == 'Iqxy': 
    8183            if not partype['orientation'] and not partype['magnetic']: 
    … …  
    9698            #self._theory = np.zeros_like(self.Iq) 
    9799            q_vectors = res.q_calc 
     100            q_mono = [] 
    98101        elif self.data_type == 'Iq': 
    99102            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    … …  
    120123            #self._theory = np.zeros_like(self.Iq) 
    121124            q_vectors = [res.q_calc] 
     125            q_mono = [] 
    122126        else: 
    123127            raise ValueError("Unknown data type") # never gets here 
    … …  
    125129        # Remember function inputs so we can delay loading the function and 
    126130        # so we can save/restore state 
    127         self._kernel_inputs = [v for v in q_vectors] 
     131        self._kernel_inputs = q_vectors 
     132        self._kernel_mono_inputs = q_mono 
    128133        self._kernel = None 
    129134        self.Iq, self.dIq, self.index = Iq, dIq, index 
    … …  
    149154    def _calc_theory(self, pars, cutoff=0.0): 
    150155        if self._kernel is None: 
    151             self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init 
     156            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
     157            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None 
    152158 
    153159        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
     160        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None 
    154161        if self.data_type == 'sesans': 
    155             result = sesans.hankel(self._data.x, self._data.lam * 1e-9, 
    156                                    self._data.sample.thickness / 10, 
    157                                    self._kernel_inputs[0], Iq_calc) 
     162            result = sesans.transform(self._data, 
     163                                   self._kernel_inputs[0], Iq_calc,  
     164                                   self._kernel_mono_inputs, Iq_mono) 
    158165        else: 
    159166            result = self.resolution.apply(Iq_calc) 
    160         return result 
     167        return result         
    161168 
    162169 

  • sasmodels/models/adsorbed_layer.py

    rf10bc52 rc079f50  
    66 
    77r""" 
    8 This model describes the scattering from a layer of surfactant or polymer adsorbed on spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for. 
     8This model describes the scattering from a layer of surfactant or polymer adsorbed on large, smooth, notionally spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for. 
    99 
    1010Unlike many other core-shell models, this model does not assume any form for the density distribution of the adsorbed species normal to the interface (cf, a core-shell model normally assumes the density distribution to be a homogeneous step-function). For comparison, if the thickness of a (traditional core-shell like) step function distribution is *t*, the second moment about the mean of the density distribution (ie, the distance of the centre-of-mass of the distribution from the interface), |sigma| =  sqrt((*t* :sup:`2` )/12). 
    … …  
    3434 
    3535description =  """ 
    36     Evaluates the scattering from particles 
     36    Evaluates the scattering from large particles 
    3737    with an adsorbed layer of surfactant or 
    3838    polymer, independent of the form of the 
    … …  
    4242 
    4343#             ["name", "units", default, [lower, upper], "type", "description"], 
    44 parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment"], 
    45               ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount"], 
    46               ["density_poly", "g/cm3", 0.7, [0.0, inf], "", "Polymer density"], 
    47               ["radius", "Ang", 500.0, [0.0, inf], "", "Particle radius"], 
    48               ["volfraction", "None", 0.14, [0.0, inf], "", "Particle vol fraction"], 
    49               ["polymer_sld", "1/Ang^2", 1.5e-06, [-inf, inf], "", "Polymer SLD"], 
    50               ["solvent_sld", "1/Ang^2", 6.3e-06, [-inf, inf], "", "Solvent SLD"]] 
     44parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment of polymer distribution"], 
     45              ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount of polymer"], 
     46              ["density_shell", "g/cm3", 0.7, [0.0, inf], "", "Bulk density of polymer in the shell"], 
     47              ["radius", "Ang", 500.0, [0.0, inf], "", "Core particle radius"], 
     48              ["volfraction", "None", 0.14, [0.0, inf], "", "Core particle volume fraction"], 
     49              ["sld_shell", "1e-6/Ang^2", 1.5, [-inf, inf], "sld", "Polymer shell SLD"], 
     50              ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent SLD"]] 
    5151 
    5252# NB: Scale and Background are implicit parameters on every model 
    53 def Iq(q, second_moment, adsorbed_amount, density_poly, radius,  
    54         volfraction, polymer_sld, solvent_sld): 
     53def Iq(q, second_moment, adsorbed_amount, density_shell, radius,  
     54        volfraction, sld_shell, sld_solvent): 
    5555    # pylint: disable = missing-docstring 
    56     deltarhosqrd =  (polymer_sld - solvent_sld) * (polymer_sld - solvent_sld) 
    57     numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount) 
    58     denominator =  (q * q) * (density_poly * density_poly) * radius 
    59     eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment)) 
    60     #scale by 10^10 for units conversion to cm^-1 
    61     inten =  1.0e+10 * deltarhosqrd * ((numerator / denominator) * eterm) 
     56#    deltarhosqrd =  (sld_shell - sld_solvent) * (sld_shell - sld_solvent) 
     57#    numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount) 
     58#    denominator =  (q * q) * (density_shell * density_shell) * radius 
     59#    eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment)) 
     60#    #scale by 10^-2 for units conversion to cm^-1 
     61#    inten =  1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm) 
     62    aa =  (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell  
     63    bb = q * second_moment 
     64    #scale by 10^-2 for units conversion to cm^-1 
     65    inten =  6.0e-02 * pi * volfraction * aa * aa * exp(-bb * bb) / radius 
    6266    return inten 
    6367Iq.vectorized =  True  # Iq accepts an array of q values 
    … …  
    7175            second_moment = 23.0, 
    7276            adsorbed_amount = 1.9, 
    73             density_poly = 0.7, 
     77            density_shell = 0.7, 
    7478            radius = 500.0, 
    7579            volfraction = 0.14, 
    76             polymer_sld = 1.5e-06, 
    77             solvent_sld = 6.3e-06, 
     80            sld_shell = 1.5, 
     81            sld_solvent = 6.3, 
    7882            background = 0.0) 
    7983 
    … …  
    8286               second_moment = 'second_moment', 
    8387               adsorbed_amount = 'ads_amount', 
    84                density_poly = 'density_poly', 
     88               density_shell = 'density_poly', 
    8589               radius = 'radius_core', 
    8690               volfraction = 'volf_cores', 
    87                polymer_sld = 'sld_poly', 
    88                solvent_sld = 'sld_solv', 
     91               sld_shell = 'sld_poly', 
     92               sld_solvent = 'sld_solv', 
    8993               background = 'background') 
    9094 
    … …  
    9296tests =  [ 
    9397    [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9,  
    94      'density_poly': 0.7, 'radius': 500.0, 'volfraction': 0.14,  
    95      'polymer_sld': 1.5e-06, 'solvent_sld': 6.3e-06, 'background': 0.0}, 
     98     'density_shell': 0.7, 'radius': 500.0, 'volfraction': 0.14,  
     99     'sld_shell': 1.5, 'sld_solvent': 6.3, 'background': 0.0}, 
    96100     [0.0106939, 0.469418], [73.741, 9.65391e-53]], 
    97101    ] 
     102# ADDED by: SMK  ON: 16Mar2016  convert from sasview, check vs SANDRA, 18Mar2016 RKH some edits & renaming 

  • sasmodels/models/core_shell_ellipsoid.c

    re7678b2 r29172aa  
    88          double equat_shell, 
    99          double polar_shell, 
    10           double core_sld, 
    11           double shell_sld, 
    12           double solvent_sld); 
     10          double sld_core, 
     11          double sld_shell, 
     12          double sld_solvent); 
    1313 
    1414 
    … …  
    1818          double equat_shell, 
    1919          double polar_shell, 
    20           double core_sld, 
    21           double shell_sld, 
    22           double solvent_sld, 
     20          double sld_core, 
     21          double sld_shell, 
     22          double sld_solvent, 
    2323          double theta, 
    2424          double phi); 
    … …  
    4040          double equat_shell, 
    4141          double polar_shell, 
    42           double core_sld, 
    43           double shell_sld, 
    44           double solvent_sld) 
     42          double sld_core, 
     43          double sld_shell, 
     44          double sld_solvent) 
    4545{ 
    4646 
    … …  
    5151    double summ = 0.0;   //initialize intergral 
    5252 
    53     const double delpc = core_sld - shell_sld;    //core - shell 
    54     const double delps = shell_sld - solvent_sld; //shell - solvent 
     53    const double delpc = sld_core - sld_shell;    //core - shell 
     54    const double delps = sld_shell - sld_solvent; //shell - solvent 
    5555 
    5656    for(int i=0;i<N_POINTS_76;i++) { 
    … …  
    8181          double equat_shell, 
    8282          double polar_shell, 
    83           double core_sld, 
    84           double shell_sld, 
    85           double solvent_sld, 
     83          double sld_core, 
     84          double sld_shell, 
     85          double sld_solvent, 
    8686          double theta, 
    8787          double phi) 
    … …  
    9696    const double cyl_y = sin(theta); 
    9797 
    98     const double sldcs = core_sld - shell_sld; 
    99     const double sldss = shell_sld- solvent_sld; 
     98    const double sldcs = sld_core - sld_shell; 
     99    const double sldss = sld_shell- sld_solvent; 
    100100 
    101101    // Compute the angle btw vector q and the 
    … …  
    124124          double equat_shell, 
    125125          double polar_shell, 
    126           double core_sld, 
    127           double shell_sld, 
    128           double solvent_sld) 
     126          double sld_core, 
     127          double sld_shell, 
     128          double sld_solvent) 
    129129{ 
    130130    double intensity = core_shell_ellipsoid_kernel(q, 
    … …  
    133133           equat_shell, 
    134134           polar_shell, 
    135            core_sld, 
    136            shell_sld, 
    137            solvent_sld); 
     135           sld_core, 
     136           sld_shell, 
     137           sld_solvent); 
    138138 
    139139    return intensity; 
    … …  
    146146          double equat_shell, 
    147147          double polar_shell, 
    148           double core_sld, 
    149           double shell_sld, 
    150           double solvent_sld, 
     148          double sld_core, 
     149          double sld_shell, 
     150          double sld_solvent, 
    151151          double theta, 
    152152          double phi) 
    … …  
    159159                       equat_shell, 
    160160                       polar_shell, 
    161                        core_sld, 
    162                        shell_sld, 
    163                        solvent_sld, 
     161                       sld_core, 
     162                       sld_shell, 
     163                       sld_solvent, 
    164164                       theta, 
    165165                       phi); 

  • sasmodels/models/core_shell_ellipsoid.py

    r5111921 r29172aa  
    11r""" 
    22This model provides the form factor, $P(q)$, for a core shell ellipsoid (below) 
    3 where the form factor is normalized by the volume of the cylinder. 
     3where the form factor is normalized by the volume of the outer [CHECK]. 
    44 
    55.. math:: 
    … …  
    77    P(q) = scale * \left<f^2\right>/V + background 
    88 
    9 where the volume $V = (4/3)\pi(r_{maj}r_{min}^2)$ and the averaging $< >$ is 
     9where the volume $V = (4/3)\pi(R_{major\_outer}R_{minor\_outer}^2)$ and the averaging $< >$ is 
    1010applied over all orientations for 1D. 
    1111 
    … …  
    2222 
    2323    P(q) = \frac{scale}{V}\int_0^1 
    24         \left|F(q,r_{min},r_{max},\alpha)\right|^2d\alpha + background 
     24        \left|F(q,r_{minor\_core},r_{major\_core},\alpha) + F(q,r_{major\_outer},r_{major\_outer},\alpha)\right|^2d\alpha + background 
    2525 
    26     \left|F(q,r_{min},r_{max},\alpha)\right|=V\Delta \rho \cdot (3j_1(u)/u) 
     26    \left|F(q,r_{minor},r_{major},\alpha)\right|=(4\pi/3)r_{major}r_{minor}^2 \Delta \rho \cdot (3j_1(u)/u) 
    2727 
    28     u = q\left[ r_{maj}^2\alpha ^2 + r_{min}^2(1-\alpha ^2)\right]^{1/2} 
     28    u = q\left[ r_{major}^2\alpha ^2 + r_{minor}^2(1-\alpha ^2)\right]^{1/2} 
    2929 
    3030where 
    … …  
    4343polar core radius, *equat_shell* = $r_{min}$ (or equatorial outer radius), 
    4444and *polar_shell* = $r_{maj}$ (or polar outer radius). 
     45 
     46Note:It is the users' responsibility to ensure that shell radii are larger than  
     47the core radii, especially if both are polydisperse, in which case the 
     48core_shell_ellipsoid_xt model may be much better. 
     49 
    4550 
    4651.. note:: 
    … …  
    7782    single particle scattering amplitude. 
    7883    [Parameters]: 
    79     equat_core = equatorial radius of core, 
    80     polar_core = polar radius of core, 
    81     equat_shell = equatorial radius of shell, 
    82     polar_shell = polar radius (revolution axis) of shell, 
    83     core_sld = SLD_core 
    84     shell_sld = SLD_shell 
    85     solvent_sld = SLD_solvent 
     84    equat_core = equatorial radius of core, Rminor_core, 
     85    polar_core = polar radius of core, Rmajor_core, 
     86    equat_shell = equatorial radius of shell, Rminor_outer, 
     87    polar_shell = polar radius of shell, Rmajor_outer, 
     88    sld_core = scattering length density of core, 
     89    sld_shell = scattering length density of shell, 
     90    sld_solvent = scattering length density of solvent, 
    8691    background = Incoherent bkg 
    8792    scale =scale 
    8893    Note:It is the users' responsibility to ensure 
    89     that shell radii are larger than core radii. 
     94    that shell radii are larger than core radii, 
     95    especially if both are polydisperse. 
    9096    oblate: polar radius < equatorial radius 
    9197    prolate :  polar radius > equatorial radius 
    … …  
    97103#             ["name", "units", default, [lower, upper], "type", "description"], 
    98104parameters = [ 
    99     ["equat_core",  "Ang",      200,   [0, inf],    "volume",      "Equatorial radius of core"], 
    100     ["polar_core",  "Ang",       10,   [0, inf],    "volume",      "Polar radius of core"], 
    101     ["equat_shell", "Ang",      250,   [0, inf],    "volume",      "Equatorial radius of shell"], 
    102     ["polar_shell", "Ang",       30,   [0, inf],    "volume",      "Polar radius of shell"], 
    103     ["core_sld",    "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"], 
    104     ["shell_sld",   "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"], 
    105     ["solvent_sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"], 
     105    ["equat_core",  "Ang",      200,   [0, inf],    "volume",      "Equatorial radius of core, Rminor_core "], 
     106    ["polar_core",  "Ang",       10,   [0, inf],    "volume",      "Polar radius of core, Rmajor_core"], 
     107    ["equat_shell", "Ang",      250,   [0, inf],    "volume",      "Equatorial radius of shell, Rminor_outer "], 
     108    ["polar_shell", "Ang",       30,   [0, inf],    "volume",      "Polar radius of shell, Rmajor_outer"], 
     109    ["sld_core",    "1e-6/Ang^2", 2,   [-inf, inf], "sld",            "Core scattering length density"], 
     110    ["sld_shell",   "1e-6/Ang^2", 1,   [-inf, inf], "sld",            "Shell scattering length density"], 
     111    ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",            "Solvent scattering length density"], 
    106112    ["theta",       "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"], 
    107113    ["phi",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"], 
    … …  
    116122            equat_shell=250.0, 
    117123            polar_shell=30.0, 
    118             core_sld=2.0, 
    119             shell_sld=1.0, 
    120             solvent_sld=6.3, 
     124            sld_core=2.0, 
     125            sld_shell=1.0, 
     126            sld_solvent=6.3, 
    121127            theta=0, 
    122128            phi=0) 
    … …  
    124130oldname = 'CoreShellEllipsoidModel' 
    125131 
    126 oldpars = dict(core_sld='sld_core', 
    127                shell_sld='sld_shell', 
    128                solvent_sld='sld_solvent', 
     132oldpars = dict(equat_core='equat_core',polar_core='polar_core',equat_shell='equat_shell',polar_shell='polar_shell', 
     133               sld_core='sld_core', 
     134               sld_shell='sld_shell', 
     135               sld_solvent='sld_solvent', 
    129136               theta='axis_theta', 
    130137               phi='axis_phi') 
    … …  
    141148      'equat_shell': 250.0, 
    142149      'polar_shell': 30.0, 
    143       'core_sld': 2.0, 
    144       'shell_sld': 1.0, 
    145       'solvent_sld': 6.3, 
     150      'sld_core': 2.0, 
     151      'sld_shell': 1.0, 
     152      'sld_solvent': 6.3, 
    146153      'background': 0.001, 
    147154      'scale': 1.0, 
    … …  
    155162      'equat_shell': 54.0, 
    156163      'polar_shell': 3.0, 
    157       'core_sld': 20.0, 
    158       'shell_sld': 10.0, 
    159       'solvent_sld': 6.0, 
     164      'sld_core': 20.0, 
     165      'sld_shell': 10.0, 
     166      'sld_solvent': 6.0, 
    160167      'background': 0.0, 
    161168      'scale': 1.0, 
    … …  
    168175      'equat_shell': 54.0, 
    169176      'polar_shell': 3.0, 
    170       'core_sld': 20.0, 
    171       'shell_sld': 10.0, 
    172       'solvent_sld': 6.0, 
     177      'sld_core': 20.0, 
     178      'sld_shell': 10.0, 
     179      'sld_solvent': 6.0, 
    173180      'background': 0.01, 
    174181      'scale': 0.01, 

  • sasmodels/models/core_shell_ellipsoid_xt.py

    r5111921 r29172aa  
    11r""" 
    2 An alternative version of $P(q)$ for the core-shell ellipsoid 
    3 (see CoreShellEllipsoidModel), having as parameters the core axial ratio X 
    4 and a shell thickness, which are more often what we would like to determine. 
     2An alternative version of $P(q)$ for the core_shell_ellipsoid 
     3having as parameters the core axial ratio X and a shell thickness,  
     4which are more often what we would like to determine. 
    55 
    66This model is also better behaved when polydispersity is applied than the four 
    7 independent radii in CoreShellEllipsoidModel. 
     7independent radii in core_shell_ellipsoid model. 
    88 
    99Definition 
    … …  
    5252---------- 
    5353 
    54 R K Heenan, *Private communication* 
     54R K Heenan, 2015, reparametrised the core_shell_ellipsoid model 
    5555 
    5656""" 
    … …  
    6969        [Parameters]: 
    7070        equat_core = equatorial radius of core, 
    71         x_core = polar radius of core, 
     71        x_core = ratio of core polar/equatorial radii, 
    7272        t_shell = equatorial radius of outer surface, 
    73         x_polar_shell = polar radius (revolution axis) of outer surface 
    74         core_sld = SLD_core 
    75         shell_sld = SLD_shell 
    76         solvent_sld = SLD_solvent 
     73        x_polar_shell = ratio of polar shell thickness to equatorial shell thickness, 
     74        sld_core = SLD_core 
     75        sld_shell = SLD_shell 
     76        sld_solvent = SLD_solvent 
    7777        background = Incoherent bkg 
    7878        scale =scale 
    … …  
    9292    ["t_shell",       "Ang",       30,   [0, inf],    "volume",      "Equatorial radius of shell"], 
    9393    ["x_polar_shell", "",           1,   [0, inf],    "volume",      "Polar radius of shell"], 
    94     ["core_sld",      "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"], 
    95     ["shell_sld",     "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"], 
    96     ["solvent_sld",   "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"], 
     94    ["sld_core",      "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"], 
     95    ["sld_shell",     "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"], 
     96    ["sld_solvent",   "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"], 
    9797    ["theta",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"], 
    9898    ["phi",           "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"], 
    … …  
    108108            t_shell=30.0, 
    109109            x_polar_shell=1.0, 
    110             core_sld=2.0, 
    111             shell_sld=1.0, 
    112             solvent_sld=6.3, 
     110            sld_core=2.0, 
     111            sld_shell=1.0, 
     112            sld_solvent=6.3, 
    113113            theta=0, 
    114114            phi=0) 
    … …  
    119119               x_core='X_core', 
    120120               t_shell='T_shell', 
    121                core_sld='sld_core', 
    122                shell_sld='sld_shell', 
    123                solvent_sld='sld_solvent', 
     121               sld_core='sld_core', 
     122               sld_shell='sld_shell', 
     123               sld_solvent='sld_solvent', 
    124124               theta='axis_theta', 
    125125               phi='axis_phi') 
    … …  
    136136      't_shell': 50.0, 
    137137      'x_polar_shell': 0.2, 
    138       'core_sld': 2.0, 
    139       'shell_sld': 1.0, 
    140       'solvent_sld': 6.3, 
     138      'sld_core': 2.0, 
     139      'sld_shell': 1.0, 
     140      'sld_solvent': 6.3, 
    141141      'background': 0.001, 
    142142      'scale': 1.0, 
    … …  
    150150      't_shell': 54.0, 
    151151      'x_polar_shell': 3.0, 
    152       'core_sld': 20.0, 
    153       'shell_sld': 10.0, 
    154       'solvent_sld': 6.0, 
     152      'sld_core': 20.0, 
     153      'sld_shell': 10.0, 
     154      'sld_solvent': 6.0, 
    155155      'background': 0.0, 
    156156      'scale': 1.0, 
    … …  
    163163      't_shell': 54.0, 
    164164      'x_polar_shell': 3.0, 
    165       'core_sld': 20.0, 
    166       'shell_sld': 10.0, 
    167       'solvent_sld': 6.0, 
     165      'sld_core': 20.0, 
     166      'sld_shell': 10.0, 
     167      'sld_solvent': 6.0, 
    168168      'background': 0.01, 
    169169      'scale': 0.01, 

  • sasmodels/models/hardsphere.py

    r3bcd03d re98c1e0  
    5858category = "structure-factor" 
    5959structure_factor = True 
     60single = False 
    6061 
    6162#             ["name", "units", default, [lower, upper], "type","description"], 

  • sasmodels/models/lamellar.py

    raa2edb2 r7c391dd  
    55---------- 
    66 
    7 The scattering intensity $I(q)$ is 
     7The scattering intensity $I(q)$ for dilute, randomly oriented, "infinitely large" sheets or lamellae is 
    88 
    99.. math:: 
    1010 
    11     I(q) = \frac{2\pi P(q)}{\delta q^2} 
     11    I(q) = scale*\frac{2\pi P(q)}{q^2\delta } 
    1212 
    1313 
    … …  
    1616.. math:: 
    1717 
    18     P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta)) 
     18   P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta)) = \frac{4\Delta\rho^2}{q^2}sin^2(\frac{q\delta}{2}) 
    1919 
     20where $\delta$ is the total layer thickness and $\Delta\rho$ is the scattering length density difference. 
    2021 
    21 where $\delta$ is the bilayer thickness. 
     22This is the limiting form for a spherical shell of infinitely large radius. Note that the division by $\delta$ 
     23means that $scale$ in sasview is the volume fraction of sheet, $\phi = S\delta$ where $S$ is the area of  
     24sheet per unit volume. $S$ is half the Porod surface area per unit volume of a thicker layer (as that would  
     25include both faces of the sheet). 
    2226 
    2327The 2D scattering intensity is calculated in the same way as 1D, where 
    … …  
    5660 
    5761#             ["name", "units", default, [lower, upper], "type","description"], 
    58 parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "", 
    59                "Layer scattering length density" ], 
    60               ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "", 
    61                "Solvent scattering length density" ], 
    62               ["thickness", "Ang", 50, [0, inf], "volume","Bilayer thickness" ], 
     62parameters = [ ["thickness", "Ang", 50, [0, inf], "volume","total layer thickness" ], 
     63               ["sld", "1e-6/Ang^2", 1, [-inf, inf], "","Layer scattering length density" ], 
     64               ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "","Solvent scattering length density" ], 
    6365             ] 
    6466 
    65  
    6667# No volume normalization despite having a volume parameter 
    67 # This should perhaps be volume normalized? 
     68# This should perhaps be volume normalized? - it is! 
    6869form_volume = """ 
    6970    return 1.0; 
    … …  
    7172 
    7273Iq = """ 
    73     const double sub = sld - solvent_sld; 
     74    const double sub = sld - sld_solvent; 
    7475    const double qsq = q*q; 
    7576    // Original expression 
    … …  
    8990 
    9091demo = dict(scale=1, background=0, 
    91             sld=6, solvent_sld=1, 
     92            sld=6, sld_solvent=1, 
    9293            thickness=40, 
    9394            thickness_pd=0.2, thickness_pd_n=40) 
    9495oldname = 'LamellarModel' 
    95 oldpars = dict(sld='sld_bi', solvent_sld='sld_sol', thickness='bi_thick') 
     96oldpars = dict(sld='sld_bi', sld_solvent='sld_sol', thickness='bi_thick') 
    9697tests = [ 
    97         [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 50.0, 'sld' : 1.0,'solvent_sld' : 6.3, 'thickness_pd' : 0.0,  
     98        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 50.0, 'sld' : 1.0,'sld_solvent' : 6.3, 'thickness_pd' : 0.0,  
    9899           }, [0.001], [882289.54309]] 
    99100        ] 

  • sasmodels/models/lamellar_hg.py

    raa2edb2 r52a3db3  
    2626 
    2727where $\delta_T$ is *tail_length*, $\delta_H$ is *head_length*, 
    28 $\Delta\rho_H$ is the head contrast (*head_sld* $-$ *solvent_sld*), 
    29 and $\Delta\rho_T$ is tail contrast (*sld* $-$ *solvent_sld*). 
     28$\Delta\rho_H$ is the head contrast (*sld_head* $-$ *sld_solvent*), 
     29and $\Delta\rho_T$ is tail contrast (*sld* $-$ *sld_solvent*). 
     30 
     31The total thickness of the lamellar sheet is $\delta_H$ + $\delta_T$ + $\delta_T$ + $\delta_H$. 
     32Note that in a non aqueous solvent the chemical "head" group may be the  
     33"Tail region" and vice-versa. 
    3034 
    3135The 2D scattering intensity is calculated in the same way as 1D, where 
    … …  
    4953from numpy import inf 
    5054 
    51 name = "lamellar_FFHG" 
    52 title = "Random lamellar phase with Head Groups" 
     55name = "lamellar_hg" 
     56title = "Random lamellar phase with Head and Tail Groups" 
    5357description = """\ 
    54     [Random lamellar phase with Head Groups] 
     58    [Random lamellar phase with Head and Tail Groups] 
    5559        I(q)= 2*pi*P(q)/(2(H+T)*q^(2)), where 
    5660        P(q)= see manual 
    5761        layer thickness =(H+T+T+H) = 2(Head+Tail) 
    5862        sld = Tail scattering length density 
    59         head_sld = Head scattering length density 
    60         solvent_sld = solvent scattering length density 
     63        sld_head = Head scattering length density 
     64        sld_solvent = solvent scattering length density 
    6165        background = incoherent background 
    6266        scale = scale factor 
    … …  
    6670# pylint: disable=bad-whitespace, line-too-long 
    6771#             ["name", "units", default, [lower, upper], "type","description"], 
    68 parameters = [["tail_length", "Ang",       15,   [0, inf],  "volume",  "Tail thickness"], 
    69               ["head_length", "Ang",       10,   [0, inf],  "volume",  "head thickness"], 
     72parameters = [["tail_length", "Ang",       15,   [0, inf],  "volume",  "Tail thickness ( total = H+T+T+H)"], 
     73              ["head_length", "Ang",       10,   [0, inf],  "volume",  "Head thickness"], 
    7074              ["sld",         "1e-6/Ang^2", 0.4, [-inf,inf], "",       "Tail scattering length density"], 
    71               ["head_sld",    "1e-6/Ang^2", 3.0, [-inf,inf], "",       "Head scattering length density"], 
    72               ["solvent_sld", "1e-6/Ang^2", 6,   [-inf,inf], "",       "Solvent scattering length density"]] 
     75              ["sld_head",    "1e-6/Ang^2", 3.0, [-inf,inf], "",       "Head scattering length density"], 
     76              ["sld_solvent", "1e-6/Ang^2", 6,   [-inf,inf], "",       "Solvent scattering length density"]] 
    7377# pylint: enable=bad-whitespace, line-too-long 
    7478 
    … …  
    8185Iq = """ 
    8286    const double qsq = q*q; 
    83     const double drh = head_sld - solvent_sld; 
    84     const double drt = sld - solvent_sld;    //correction 13FEB06 by L.Porcar 
     87    const double drh = sld_head - sld_solvent; 
     88    const double drt = sld - sld_solvent;    //correction 13FEB06 by L.Porcar 
    8589    const double qT = q*tail_length; 
    8690    double Pq, inten; 
    … …  
    106110demo = dict(scale=1, background=0, 
    107111            tail_length=15, head_length=10, 
    108             sld=0.4, head_sld=3.0, solvent_sld=6.0, 
     112            sld=0.4, sld_head=3.0, sld_solvent=6.0, 
    109113            tail_length_pd=0.2, tail_length_pd_n=40, 
    110114            head_length_pd=0.01, head_length_pd_n=40) 
    … …  
    112116oldname = 'LamellarFFHGModel' 
    113117oldpars = dict(head_length='h_thickness', tail_length='t_length', 
    114                sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent') 
     118               sld='sld_tail', sld_head='sld_head', sld_solvent='sld_solvent') 
    115119# 
    116120tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 15.0, 'head_length': 10.0, 
    117            'sld': 0.4, 'head_sld': 3.0, 'solvent_sld': 6.0}, [0.001], [653143.9209]]] 
     121           'sld': 0.4, 'sld_head': 3.0, 'sld_solvent': 6.0}, [0.001], [653143.9209]]] 
     122# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs 

  • sasmodels/models/lamellar_hg_stack_caille.py

    raa2edb2 r6ab4ed8  
    99.. math:: 
    1010 
    11     I(q) = 2 \pi \frac{P(q)S(q)}{\delta q^2} 
     11    I(q) = 2 \pi \frac{P(q)S(q)}{q^2\delta } 
    1212 
    1313 
    … …  
    5656results for the next lower and higher values. 
    5757 
     58Be aware that the computations may be very slow. 
     59 
    5860The 2D scattering intensity is calculated in the same way as 1D, where 
    5961the $q$ vector is defined as 
    … …  
    7375from numpy import inf 
    7476 
    75 name = "lamellarCailleHG" 
    76 title = "Random lamellar sheet with Caille structure factor" 
     77name = "lamellar_hg_stack_caille" 
     78title = "Random lamellar head/tail/tail/head sheet with Caille structure factor" 
    7779description = """\ 
    7880    [Random lamellar phase with Caille  structure factor] 
    … …  
    104106    ["sld", "1e-6/Ang^2", 0.4, [-inf, inf], "", 
    105107     "Tail scattering length density"], 
    106     ["head_sld", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
     108    ["sld_head", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
    107109     "Head scattering length density"], 
    108     ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "", 
     110    ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "", 
    109111     "Solvent scattering length density"], 
    110112    ] 
    111113 
    112 source = ["lamellarCailleHG_kernel.c"] 
     114source = ["lamellar_hg_stack_caille_kernel.c"] 
    113115 
    114116# No volume normalization despite having a volume parameter 
    … …  
    129131    Nlayers=20, spacing=200., Caille_parameter=0.05, 
    130132    tail_length=15, head_length=10, 
    131     #sld=-1, head_sld=4.0, solvent_sld=6.0, 
    132     sld=-1, head_sld=4.1, solvent_sld=6.0, 
     133    #sld=-1, sld_head=4.0, sld_solvent=6.0, 
     134    sld=-1, sld_head=4.1, sld_solvent=6.0, 
    133135    tail_length_pd=0.1, tail_length_pd_n=20, 
    134136    head_length_pd=0.05, head_length_pd_n=30, 
    … …  
    140142oldpars = dict( 
    141143    tail_length='deltaT', head_length='deltaH', Nlayers='n_plates', 
    142     Caille_parameter='caille', sld='sld_tail', head_sld='sld_head', 
    143     solvent_sld='sld_solvent') 
     144    Caille_parameter='caille', sld='sld_tail', sld_head='sld_head', 
     145    sld_solvent='sld_solvent') 
    144146# 
    145147tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 10.0, 'head_length': 2.0, 
    146148           'Nlayers': 30.0, 'spacing': 40., 'Caille_parameter': 0.001, 'sld': 0.4, 
    147            'head_sld': 2.0, 'solvent_sld': 6.0, 'tail_length_pd': 0.0, 
     149           'sld_head': 2.0, 'sld_solvent': 6.0, 'tail_length_pd': 0.0, 
    148150           'head_length_pd': 0.0, 'spacing_pd': 0.0}, [0.001], [6838238.571488]]] 
     151# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs 

  • sasmodels/models/lamellar_stack_caille.py

    raa2edb2 r6ab4ed8  
    1111.. math:: 
    1212 
    13     I(q) = 2\pi \frac{P(q)S(q)}{\delta q^2} 
     13    I(q) = 2\pi \frac{P(q)S(q)}{q^2\delta } 
    1414 
    1515The form factor is 
    … …  
    7070from numpy import inf 
    7171 
    72 name = "lamellarPS" 
     72name = "lamellar_stack_caille" 
    7373title = "Random lamellar sheet with Caille structure factor" 
    7474description = """\ 
    … …  
    9292    ["Caille_parameter", "1/Ang^2",    0.1, [0.0,0.8],  "",       "Caille parameter"], 
    9393    ["sld",              "1e-6/Ang^2", 6.3, [-inf,inf], "",       "layer scattering length density"], 
    94     ["solvent_sld",      "1e-6/Ang^2", 1.0, [-inf,inf], "",       "Solvent scattering length density"], 
     94    ["sld_solvent",      "1e-6/Ang^2", 1.0, [-inf,inf], "",       "Solvent scattering length density"], 
    9595    ] 
    9696# pylint: enable=bad-whitespace, line-too-long 
    9797 
    98 source = ["lamellarCaille_kernel.c"] 
     98source = ["lamellar_stack_caille_kernel.c"] 
    9999 
    100100# No volume normalization despite having a volume parameter 
    … …  
    113113demo = dict(scale=1, background=0, 
    114114            thickness=67., Nlayers=3.75, spacing=200., 
    115             Caille_parameter=0.268, sld=1.0, solvent_sld=6.34, 
     115            Caille_parameter=0.268, sld=1.0, sld_solvent=6.34, 
    116116            thickness_pd=0.1, thickness_pd_n=100, 
    117117            spacing_pd=0.05, spacing_pd_n=40) 
    … …  
    120120oldpars = dict(thickness='delta', Nlayers='N_plates', 
    121121               Caille_parameter='caille', 
    122                sld='sld_bi', solvent_sld='sld_sol') 
     122               sld='sld_bi', sld_solvent='sld_sol') 
    123123# 
    124124tests = [ 
    125125    [{'scale': 1.0, 'background': 0.0, 'thickness': 30., 'Nlayers': 20.0, 
    126126      'spacing': 400., 'Caille_parameter': 0.1, 'sld': 6.3, 
    127       'solvent_sld': 1.0, 'thickness_pd': 0.0, 'spacing_pd': 0.0}, 
     127      'sld_solvent': 1.0, 'thickness_pd': 0.0, 'spacing_pd': 0.0}, 
    128128     [0.001], [28895.13397]] 
    129129    ] 
     130# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs 

  • sasmodels/models/lamellar_stack_paracrystal.py

    raa2edb2 r6ab4ed8  
    1616  *not* the total excluded volume of the paracrystal), 
    1717 
    18 - *sld* $-$ *solvent_sld* is the contrast $\Delta \rho$, 
     18- *sld* $-$ *sld_solvent* is the contrast $\Delta \rho$, 
    1919 
    2020- *thickness* is the layer thickness $t$, 
    … …  
    3636 
    3737The form factor of the bilayer is approximated as the cross section of an 
    38 infinite, planar bilayer of thickness $t$ 
     38infinite, planar bilayer of thickness $t$ (compare the equations for the 
     39lamellar model). 
    3940 
    4041.. math:: 
    … …  
    9495from numpy import inf 
    9596 
    96 name = "lamellarPC" 
     97name = "lamellar_stack_paracrystal" 
    9798title = "Random lamellar sheet with paracrystal structure factor" 
    9899description = """\ 
    … …  
    120121              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
    121122               "layer scattering length density"], 
    122               ["solvent_sld", "1e-6/Ang^2", 6.34, [-inf, inf], "", 
     123              ["sld_solvent", "1e-6/Ang^2", 6.34, [-inf, inf], "", 
    123124               "Solvent scattering length density"], 
    124125             ] 
    125126 
    126127 
    127 source = ["lamellarPC_kernel.c"] 
     128source = ["lamellar_stack_paracrystal_kernel.c"] 
    128129 
    129130form_volume = """ 
    … …  
    140141demo = dict(scale=1, background=0, 
    141142            thickness=33, Nlayers=20, spacing=250, spacing_polydisp=0.2, 
    142             sld=1.0, solvent_sld=6.34, 
     143            sld=1.0, sld_solvent=6.34, 
    143144            thickness_pd=0.2, thickness_pd_n=40) 
    144145 
    145146oldname = 'LamellarPCrystalModel' 
    146147oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer', 
    147                solvent_sld='sld_solvent') 
     148               sld_solvent='sld_solvent') 
    148149# 
    149150tests = [ 
    150151        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 33.,'Nlayers' : 20.0, 'spacing' : 250., 'spacing_polydisp' : 0.2, 
    151             'sld' : 1.0, 'solvent_sld' : 6.34, 
     152            'sld' : 1.0, 'sld_solvent' : 6.34, 
    152153            'thickness_pd' : 0.0, 'thickness_pd_n' : 40 }, [0.001, 0.215268], [21829.3, 0.00487686]] 
    153154        ] 
     155# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs 

  • sasmodels/product.py

    r3bcd03d rd5ba841  
    124124        # so borrow values from end of p_fixed.  This makes volfraction the 
    125125        # first S parameter. 
    126         start += num_p_fixed - 2 
    127         par_map['s_fixed'] = np.arange(start, start+num_s_fixed) 
     126        start += num_p_fixed 
     127        par_map['s_fixed'] = np.hstack(([start,start], 
     128                                        np.arange(start, start+num_s_fixed-2))) 
    128129        par_map['volfraction'] = num_p_fixed 
    129         start += num_s_fixed 
     130        start += num_s_fixed-2 
    130131        # vol pars offset from the start of pd pars 
    131132        par_map['vol_pars'] = [start+k for k in vol_par_idx] 
    132133        par_map['p_pd'] = np.arange(start, start+num_p_pd) 
    133         start += num_p_pd 
    134         par_map['s_pd'] = np.arange(start-1, start+num_s_pd)  # should be empty... 
     134        start += num_p_pd-1 
     135        par_map['s_pd'] = np.hstack((start, 
     136                                     np.arange(start, start+num_s_pd-1))) 
    135137 
    136138        self.fixed_pars = model_info['partype']['fixed-' + dim] 

  • sasmodels/resolution.py

    r17bbadd ra146eaa  
    197197 
    198198 
    199     **Algorithm** 
     199    Definition 
     200    ---------- 
    200201 
    201202    We are using the mid-point integration rule to assign weights to each 
    … …  
    441442    .. math:: 
    442443 
    443         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n) 
     444         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n) 
    444445            / (\log q_n - log q_1) 
    445446    """ 

  • sasmodels/sasview_model.py

    rfcd7bbd r2622b3f  
    1717from copy import deepcopy 
    1818import warnings 
     19import collections 
    1920 
    2021import numpy as np 
    … …  
    5758        ## TODO: reorganize parameter handling 
    5859        self.details = dict() 
    59         self.params = dict() 
     60        self.params = collections.OrderedDict() 
    6061        self.dispersion = dict() 
    6162        partype = model.info['partype'] 
     63 
    6264        for p in model.info['parameters']: 
    6365            self.params[p.name] = p.default 

  • sasmodels/sesans.py

    r190fc2b ra06430c  
    1313from numpy import pi, exp 
    1414from scipy.special import jv as besselj 
    15  
     15#import direct_model.DataMixin as model 
     16         
    1617def make_q(q_max, Rmax): 
    1718    r""" 
    … …  
    2122    q_min = dq = 0.1 * 2*pi / Rmax 
    2223    return np.arange(q_min, q_max, dq) 
     24     
     25def make_all_q(data): 
     26    if not data.needs_all_q: 
     27        return [] 
     28    elif needs_Iqxy(data): 
     29        # compute qx, qy 
     30        Qx, Qy = np.meshgrid(qx, qy) 
     31        return [Qx, Qy] 
     32    else: 
     33        # else only need q 
     34        return [q] 
    2335 
     36def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
     37    nqmono = len(qmono) 
     38    if nqmono == 0: 
     39        result = call_hankel(data, q_calc, Iq_calc) 
     40    elif nqmono == 1: 
     41        q = qmono[0] 
     42        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
     43    else: 
     44        Qx, Qy = [qmono[0], qmono[1]] 
     45        Qx = np.reshape(Qx, nqx, nqy) 
     46        Qy = np.reshape(Qy, nqx, nqy) 
     47        Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
     48        qx = Qx[0, :] 
     49        qy = Qy[:, 0] 
     50        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     51 
     52    return result 
     53 
     54def call_hankel(data, q_calc, Iq_calc): 
     55    return hankel(data.x, data.lam * 1e-9, 
     56                  data.sample.thickness / 10, 
     57                  q_calc, Iq_calc) 
     58   
     59def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
     60    return hankel(data.x, data.lam * 1e-9, 
     61                  data.sample.thickness / 10, 
     62                  q_calc, Iq_calc) 
     63                   
     64def Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
     65    return hankel(data.x, data.y, data.lam * 1e-9, 
     66                  data.sample.thickness / 10, 
     67                  q_calc, Iq_calc) 
     68                         
     69def TotalScatter(model, parameters):  #Work in progress!! 
     70#    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
     71    allq = np.linspace(0,4*pi/wavelength,1000) 
     72    allIq = 1 
     73    integral = allq*allIq 
     74     
     75 
     76 
     77def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
     78#============================================================================== 
     79#     2D Cosine Transform if "wavelength" is a vector 
     80#============================================================================== 
     81#allq is the q-space needed to create the total scattering cross-section 
     82 
     83    Gprime = np.zeros_like(wavelength, 'd') 
     84    s = np.zeros_like(wavelength, 'd') 
     85    sd = np.zeros_like(wavelength, 'd') 
     86    Gprime = np.zeros_like(wavelength, 'd') 
     87    f = np.zeros_like(wavelength, 'd') 
     88    for i, wavelength_i in enumerate(wavelength): 
     89        z = magfield*wavelength_i 
     90        allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
     91        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
     92        alldq = (allq[1]-allq[0])*1e10 
     93        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
     94        s[i]=1-exp(-sigma) 
     95        for j, Iqy_j, qy_j in enumerate(qy): 
     96            for k, Iqz_k, qz_k in enumerate(qz): 
     97                Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
     98                q = np.sqrt(qy_j^2 + qz_k^2) 
     99                Gintegral = Iq*cos(z*Qz_k) 
     100                Gprime[i] += Gintegral 
     101#                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
     102#                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
     103#                For now, work with standard 2-phase scatter 
     104 
     105 
     106                sd[i] += Iq 
     107        f[i] = 1-s[i]+sd[i] 
     108        P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
     109 
     110 
     111 
     112 
     113def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
     114#============================================================================== 
     115#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
     116#============================================================================== 
     117#acceptq is the q-space needed to create limited acceptance effect 
     118    SElength= wavelength*magfield 
     119    G = np.zeros_like(SElength, 'd') 
     120    threshold=2*pi*theta/wavelength 
     121    for i, SElength_i in enumerate(SElength): 
     122        allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
     123        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
     124        alldq = (allq[1]-allq[0])*1e10 
     125        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
     126        s[i]=1-exp(-sigma) 
     127 
     128        dq = (q[1]-q[0])*1e10 
     129        a = (x<threshold) 
     130        acceptq = a*q 
     131        acceptIq = a*Iq 
     132 
     133        G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
     134 
     135#        G[i]=np.sum(integral) 
     136 
     137    G *= dq*1e10*2*pi 
     138 
     139    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
     140     
    24141def hankel(SElength, wavelength, thickness, q, Iq): 
    25142    r""" 
    … …  
    44161    """ 
    45162    G = np.zeros_like(SElength, 'd') 
     163#============================================================================== 
     164#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
     165#============================================================================== 
    46166    for i, SElength_i in enumerate(SElength): 
    47167        integral = besselj(0, q*SElength_i)*Iq*q 

  • setup.py

    r3eb3312 rf903f0a  
    33packages = find_packages(exclude=['contrib', 'docs', 'tests*']) 
    44package_data = { 
    5     'sasmodels.models': ['*.c'], 
     5    'sasmodels.models': ['*.c','lib/*.c'], 
     6    'sasmodels': ['*.c'], 
    67} 
    78required = [] 

  • sasmodels/models/bessel.py

    rcbd37a7 ra5af4e1  
    7878 
    7979Iq = """ 
    80     return 2.0*j1(q)/q; 
     80    return J1(q); 
    8181    """ 
    8282 

  • sasmodels/models/lib/j1_cephes.c

    re2af2a9 ra5af4e1  
    3939Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    4040*/ 
    41 double j1(double ); 
     41double J1(double ); 
    4242 
    43 double j1(double x) { 
     43double J1(double x) { 
    4444 
    4545//Cephes double pression function 
    … …  
    4848    double w, z, p, q, xn; 
    4949 
    50     const double DR1 = 5.78318596294678452118E0; 
    51     const double DR2 = 3.04712623436620863991E1; 
    5250    const double Z1 = 1.46819706421238932572E1; 
    5351    const double Z2 = 4.92184563216946036703E1; 
    5452    const double THPIO4 =  2.35619449019234492885; 
    5553    const double SQ2OPI = 0.79788456080286535588; 
    56  
    57     double RP[8] = { 
    58     -8.99971225705559398224E8, 
    59     4.52228297998194034323E11, 
    60     -7.27494245221818276015E13, 
    61     3.68295732863852883286E15, 
    62     0.0, 
    63     0.0, 
    64     0.0, 
    65     0.0 
    66     }; 
    67  
    68     double RQ[8] = { 
    69     /* 1.00000000000000000000E0,*/ 
    70     6.20836478118054335476E2, 
    71     2.56987256757748830383E5, 
    72     8.35146791431949253037E7, 
    73     2.21511595479792499675E10, 
    74     4.74914122079991414898E12, 
    75     7.84369607876235854894E14, 
    76     8.95222336184627338078E16, 
    77     5.32278620332680085395E18, 
    78     }; 
    79  
    80     double PP[8] = { 
    81     7.62125616208173112003E-4, 
    82     7.31397056940917570436E-2, 
    83     1.12719608129684925192E0, 
    84     5.11207951146807644818E0, 
    85     8.42404590141772420927E0, 
    86     5.21451598682361504063E0, 
    87     1.00000000000000000254E0, 
    88     0.0, 
    89     }; 
    90     double PQ[8] = { 
    91     5.71323128072548699714E-4, 
    92     6.88455908754495404082E-2, 
    93     1.10514232634061696926E0, 
    94     5.07386386128601488557E0, 
    95     8.39985554327604159757E0, 
    96     5.20982848682361821619E0, 
    97     9.99999999999999997461E-1, 
    98     0.0, 
    99     }; 
    100  
    101     double QP[8] = { 
    102     5.10862594750176621635E-2, 
    103     4.98213872951233449420E0, 
    104     7.58238284132545283818E1, 
    105     3.66779609360150777800E2, 
    106     7.10856304998926107277E2, 
    107     5.97489612400613639965E2, 
    108     2.11688757100572135698E2, 
    109     2.52070205858023719784E1, 
    110     }; 
    111  
    112     double QQ[8] = { 
    113     /* 1.00000000000000000000E0,*/ 
    114     7.42373277035675149943E1, 
    115     1.05644886038262816351E3, 
    116     4.98641058337653607651E3, 
    117     9.56231892404756170795E3, 
    118     7.99704160447350683650E3, 
    119     2.82619278517639096600E3, 
    120     3.36093607810698293419E2, 
    121     0.0, 
    122     }; 
    12354 
    12455    w = x; 
    … …  
    12960        { 
    13061            z = x * x; 
    131             w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 
     62            w = polevlRP( z, 3 ) / p1evlRQ( z, 8 ); 
    13263            w = w * x * (z - Z1) * (z - Z2); 
    13364            return( w ); 
    … …  
    13667    w = 5.0/x; 
    13768    z = w * w; 
    138     p = polevl( z, PP, 6)/polevl( z, PQ, 6 ); 
    139     q = polevl( z, QP, 7)/p1evl( z, QQ, 7 ); 
     69 
     70    p = polevlPP( z, 6)/polevlPQ( z, 6 ); 
     71    q = polevlQP( z, 7)/p1evlQQ( z, 7 ); 
     72 
    14073    xn = x - THPIO4; 
    14174 
    … …  
    15588 
    15689 
    157     double JP[8] = { 
    158         -4.878788132172128E-009, 
    159         6.009061827883699E-007, 
    160         -4.541343896997497E-005, 
    161         1.937383947804541E-003, 
    162         -3.405537384615824E-002, 
    163         0.0, 
    164         0.0, 
    165         0.0 
    166     }; 
    167  
    168     double MO1[8] = { 
    169         6.913942741265801E-002, 
    170         -2.284801500053359E-001, 
    171         3.138238455499697E-001, 
    172         -2.102302420403875E-001, 
    173         5.435364690523026E-003, 
    174         1.493389585089498E-001, 
    175         4.976029650847191E-006, 
    176         7.978845453073848E-001 
    177     }; 
    178  
    179     double PH1[8] = { 
    180         -4.497014141919556E+001, 
    181         5.073465654089319E+001, 
    182         -2.485774108720340E+001, 
    183         7.222973196770240E+000, 
    184         -1.544842782180211E+000, 
    185         3.503787691653334E-001, 
    186         -1.637986776941202E-001, 
    187         3.749989509080821E-001 
    188     }; 
    189  
    19090    xx = x; 
    19191    if( xx < 0 ) 
    … …  
    19595        { 
    19696            z = xx * xx; 
    197             p = (z-Z1) * xx * polevl( z, JP, 4 ); 
     97            p = (z-Z1) * xx * polevlJP( z, 4 ); 
    19898            return( p ); 
    19999        } 
    … …  
    202102    w = sqrt(q); 
    203103 
    204     p = w * polevl( q, MO1, 7); 
     104    p = w * polevlMO1( q, 7); 
    205105    w = q*q; 
    206     xn = q * polevl( w, PH1, 7) - THPIO4F; 
     106    xn = q * polevlPH1( w, 7) - THPIO4F; 
    207107    p = p * cos(xn + xx); 
    208108 
    … …  
    210110#endif 
    211111} 
    212  

  • sasmodels/models/lib/polevl.c

    re2af2a9 ra5af4e1  
    5151*/ 
    5252 
    53 double polevl( double x, double coef[8], int N ); 
    54 double p1evl( double x, double coef[8], int N ); 
    55  
    56 double polevl( double x, double coef[8], int N ) { 
    57  
    58     int i = 0; 
    59     double ans = coef[i]; 
    60  
    61     while (i < N) { 
    62         i++; 
    63         ans = ans * x + coef[i]; 
    64     } 
    65  
    66     return ans ; 
    67  
    68 } 
     53double polevlRP(double x, int N ) { 
     54 
     55    double coef[8] = { 
     56    -8.99971225705559398224E8, 
     57    4.52228297998194034323E11, 
     58    -7.27494245221818276015E13, 
     59    3.68295732863852883286E15, 
     60    0.0, 
     61    0.0, 
     62    0.0, 
     63    0.0 }; 
     64 
     65    int i = 0; 
     66    double ans = coef[i]; 
     67 
     68    while (i < N) { 
     69        i++; 
     70        ans = ans * x + coef[i]; 
     71    } 
     72    return ans ; 
     73} 
     74 
     75double polevlRQ(double x, int N ) { 
     76 
     77    double coef[8] = { 
     78        6.20836478118054335476E2, 
     79        2.56987256757748830383E5, 
     80        8.35146791431949253037E7, 
     81        2.21511595479792499675E10, 
     82        4.74914122079991414898E12, 
     83        7.84369607876235854894E14, 
     84        8.95222336184627338078E16, 
     85        5.32278620332680085395E18 
     86    }; 
     87 
     88    int i = 0; 
     89    double ans = coef[i]; 
     90 
     91    while (i < N) { 
     92        i++; 
     93        ans = ans * x + coef[i]; 
     94    } 
     95    return ans ; 
     96} 
     97 
     98double polevlPP(double x, int N ) { 
     99 
     100    double coef[8] = { 
     101    7.62125616208173112003E-4, 
     102    7.31397056940917570436E-2, 
     103    1.12719608129684925192E0, 
     104    5.11207951146807644818E0, 
     105    8.42404590141772420927E0, 
     106    5.21451598682361504063E0, 
     107    1.00000000000000000254E0, 
     108    0.0} ; 
     109 
     110    int i = 0; 
     111    double ans = coef[i]; 
     112 
     113    while (i < N) { 
     114        i++; 
     115        ans = ans * x + coef[i]; 
     116    } 
     117    return ans ; 
     118} 
     119 
     120double polevlPQ(double x, int N ) { 
     121//4: PQ 
     122    double coef[8] = { 
     123    5.71323128072548699714E-4, 
     124    6.88455908754495404082E-2, 
     125    1.10514232634061696926E0, 
     126    5.07386386128601488557E0, 
     127    8.39985554327604159757E0, 
     128    5.20982848682361821619E0, 
     129    9.99999999999999997461E-1, 
     130    0.0 }; 
     131 
     132    int i = 0; 
     133    double ans = coef[i]; 
     134 
     135    while (i < N) { 
     136        i++; 
     137        ans = ans * x + coef[i]; 
     138    } 
     139    return ans ; 
     140} 
     141 
     142double polevlQP(double x, int N ) { 
     143    double coef[8] = { 
     144    5.10862594750176621635E-2, 
     145    4.98213872951233449420E0, 
     146    7.58238284132545283818E1, 
     147    3.66779609360150777800E2, 
     148    7.10856304998926107277E2, 
     149    5.97489612400613639965E2, 
     150    2.11688757100572135698E2, 
     151    2.52070205858023719784E1 }; 
     152 
     153    int i = 0; 
     154    double ans = coef[i]; 
     155 
     156    while (i < N) { 
     157        i++; 
     158        ans = ans * x + coef[i]; 
     159    } 
     160    return ans ; 
     161 
     162} 
     163 
     164double polevlQQ(double x, int N ) { 
     165    double coef[8] = { 
     166    7.42373277035675149943E1, 
     167    1.05644886038262816351E3, 
     168    4.98641058337653607651E3, 
     169    9.56231892404756170795E3, 
     170    7.99704160447350683650E3, 
     171    2.82619278517639096600E3, 
     172    3.36093607810698293419E2, 
     173    0.0 }; 
     174 
     175    int i = 0; 
     176    double ans = coef[i]; 
     177 
     178    while (i < N) { 
     179        i++; 
     180        ans = ans * x + coef[i]; 
     181    } 
     182    return ans ; 
     183 
     184} 
     185 
     186double polevlJP( double x, int N ) { 
     187    double coef[8] = { 
     188        -4.878788132172128E-009, 
     189        6.009061827883699E-007, 
     190        -4.541343896997497E-005, 
     191        1.937383947804541E-003, 
     192        -3.405537384615824E-002, 
     193        0.0, 
     194        0.0, 
     195        0.0 
     196    }; 
     197 
     198    int i = 0; 
     199    double ans = coef[i]; 
     200 
     201    while (i < N) { 
     202        i++; 
     203        ans = ans * x + coef[i]; 
     204    } 
     205    return ans ; 
     206 
     207} 
     208 
     209double polevlMO1( double x, int N ) { 
     210    double coef[8] = { 
     211        6.913942741265801E-002, 
     212        -2.284801500053359E-001, 
     213        3.138238455499697E-001, 
     214        -2.102302420403875E-001, 
     215        5.435364690523026E-003, 
     216        1.493389585089498E-001, 
     217        4.976029650847191E-006, 
     218        7.978845453073848E-001 
     219    }; 
     220 
     221    int i = 0; 
     222    double ans = coef[i]; 
     223 
     224    while (i < N) { 
     225        i++; 
     226        ans = ans * x + coef[i]; 
     227    } 
     228    return ans ; 
     229} 
     230 
     231double polevlPH1( double x, int N ) { 
     232 
     233    double coef[8] = { 
     234        -4.497014141919556E+001, 
     235        5.073465654089319E+001, 
     236        -2.485774108720340E+001, 
     237        7.222973196770240E+000, 
     238        -1.544842782180211E+000, 
     239        3.503787691653334E-001, 
     240        -1.637986776941202E-001, 
     241        3.749989509080821E-001 
     242    }; 
     243 
     244    int i = 0; 
     245    double ans = coef[i]; 
     246 
     247    while (i < N) { 
     248        i++; 
     249        ans = ans * x + coef[i]; 
     250    } 
     251    return ans ; 
     252} 
     253 
     254/*double polevl( double x, double coef[8], int N ) { 
     255 
     256    int i = 0; 
     257    double ans = coef[i]; 
     258 
     259    while (i < N) { 
     260        i++; 
     261        ans = ans * x + coef[i]; 
     262    } 
     263 
     264    return ans ; 
     265 
     266}*/ 
    69267 
    70268/*                                                      p1evl() */ 
    … …  
    74272 */ 
    75273 
    76 double p1evl( double x, double coef[8], int N ) { 
    77     int i=0; 
    78     double ans = x+coef[i]; 
    79  
    80     while (i < N-1) { 
    81         i++; 
    82         ans = ans*x + coef[i]; 
    83     } 
    84  
    85     return( ans ); 
    86  
    87 } 
     274double p1evlRP( double x, int N ) { 
     275    double coef[8] = { 
     276    -8.99971225705559398224E8, 
     277    4.52228297998194034323E11, 
     278    -7.27494245221818276015E13, 
     279    3.68295732863852883286E15, 
     280    0.0, 
     281    0.0, 
     282    0.0, 
     283    0.0 }; 
     284 
     285    int i=0; 
     286    double ans = x+coef[i]; 
     287 
     288    while (i < N-1) { 
     289        i++; 
     290        ans = ans*x + coef[i]; 
     291    } 
     292 
     293    return( ans ); 
     294 
     295} 
     296 
     297double p1evlRQ( double x, int N ) { 
     298//1: RQ 
     299    double coef[8] = { 
     300    6.20836478118054335476E2, 
     301    2.56987256757748830383E5, 
     302    8.35146791431949253037E7, 
     303    2.21511595479792499675E10, 
     304    4.74914122079991414898E12, 
     305    7.84369607876235854894E14, 
     306    8.95222336184627338078E16, 
     307    5.32278620332680085395E18}; 
     308 
     309    int i=0; 
     310    double ans = x+coef[i]; 
     311 
     312    while (i < N-1) { 
     313        i++; 
     314        ans = ans*x + coef[i]; 
     315    } 
     316 
     317    return( ans ); 
     318} 
     319 
     320double p1evlPP( double x, int N ) { 
     321//3 : PP 
     322    double coef[8] = { 
     323    7.62125616208173112003E-4, 
     324    7.31397056940917570436E-2, 
     325    1.12719608129684925192E0, 
     326    5.11207951146807644818E0, 
     327    8.42404590141772420927E0, 
     328    5.21451598682361504063E0, 
     329    1.00000000000000000254E0, 
     330    0.0}; 
     331 
     332    int i=0; 
     333    double ans = x+coef[i]; 
     334 
     335    while (i < N-1) { 
     336        i++; 
     337        ans = ans*x + coef[i]; 
     338    } 
     339 
     340    return( ans ); 
     341} 
     342 
     343double p1evlPQ( double x, int N ) { 
     344//4: PQ 
     345    double coef[8] = { 
     346    5.71323128072548699714E-4, 
     347    6.88455908754495404082E-2, 
     348    1.10514232634061696926E0, 
     349    5.07386386128601488557E0, 
     350    8.39985554327604159757E0, 
     351    5.20982848682361821619E0, 
     352    9.99999999999999997461E-1, 
     353    0.0}; 
     354 
     355    int i=0; 
     356    double ans = x+coef[i]; 
     357 
     358    while (i < N-1) { 
     359        i++; 
     360        ans = ans*x + coef[i]; 
     361    } 
     362 
     363    return( ans ); 
     364} 
     365 
     366double p1evlQP( double x, int N ) { 
     367//5: QP 
     368    double coef[8] = { 
     369    5.10862594750176621635E-2, 
     370    4.98213872951233449420E0, 
     371    7.58238284132545283818E1, 
     372    3.66779609360150777800E2, 
     373    7.10856304998926107277E2, 
     374    5.97489612400613639965E2, 
     375    2.11688757100572135698E2, 
     376    2.52070205858023719784E1 }; 
     377 
     378    int i=0; 
     379    double ans = x+coef[i]; 
     380 
     381    while (i < N-1) { 
     382        i++; 
     383        ans = ans*x + coef[i]; 
     384    } 
     385 
     386    return( ans ); 
     387} 
     388 
     389double p1evlQQ( double x, int N ) { 
     390    double coef[8] = { 
     391    7.42373277035675149943E1, 
     392    1.05644886038262816351E3, 
     393    4.98641058337653607651E3, 
     394    9.56231892404756170795E3, 
     395    7.99704160447350683650E3, 
     396    2.82619278517639096600E3, 
     397    3.36093607810698293419E2, 
     398    0.0}; 
     399 
     400    int i=0; 
     401    double ans = x+coef[i]; 
     402 
     403    while (i < N-1) { 
     404        i++; 
     405        ans = ans*x + coef[i]; 
     406    } 
     407 
     408    return( ans ); 
     409 
     410} 
     411 
     412double p1evlJP( double x, int N ) { 
     413    double coef[8] = { 
     414        -4.878788132172128E-009, 
     415        6.009061827883699E-007, 
     416        -4.541343896997497E-005, 
     417        1.937383947804541E-003, 
     418        -3.405537384615824E-002, 
     419        0.0, 
     420        0.0, 
     421        0.0}; 
     422 
     423    int i=0; 
     424    double ans = x+coef[i]; 
     425 
     426    while (i < N-1) { 
     427        i++; 
     428        ans = ans*x + coef[i]; 
     429    } 
     430 
     431    return( ans ); 
     432} 
     433 
     434double p1evlMO1( double x, int N ) { 
     435    double coef[8] = { 
     436        6.913942741265801E-002, 
     437        -2.284801500053359E-001, 
     438        3.138238455499697E-001, 
     439        -2.102302420403875E-001, 
     440        5.435364690523026E-003, 
     441        1.493389585089498E-001, 
     442        4.976029650847191E-006, 
     443        7.978845453073848E-001}; 
     444 
     445    int i=0; 
     446    double ans = x+coef[i]; 
     447 
     448    while (i < N-1) { 
     449        i++; 
     450        ans = ans*x + coef[i]; 
     451    } 
     452 
     453    return( ans ); 
     454 
     455} 
     456 
     457double p1evlPH1( double x, int N ) { 
     458    double coef[8] = { 
     459        -4.497014141919556E+001, 
     460        5.073465654089319E+001, 
     461        -2.485774108720340E+001, 
     462        7.222973196770240E+000, 
     463        -1.544842782180211E+000, 
     464        3.503787691653334E-001, 
     465        -1.637986776941202E-001, 
     466        3.749989509080821E-001}; 
     467    int i=0; 
     468    double ans = x+coef[i]; 
     469 
     470    while (i < N-1) { 
     471        i++; 
     472        ans = ans*x + coef[i]; 
     473    } 
     474 
     475    return( ans ); 
     476} 
     477 
     478/*double p1evl( double x, double coef[8], int N ) { 
     479    int i=0; 
     480    double ans = x+coef[i]; 
     481 
     482    while (i < N-1) { 
     483        i++; 
     484        ans = ans*x + coef[i]; 
     485    } 
     486 
     487    return( ans ); 
     488 
     489}*/ 
Note: See TracChangeset for help on using the changeset viewer.