Changeset 49284e1 in sasmodels for sasmodels

Timestamp:

Jan 25, 2018 12:05:09 PM (6 years ago)

Author:

GitHub <noreply@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

a69d8cd, 95498a3

Parents:

3fb9404 (diff), 7a516d0 (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.

git-author:

Adam Washington <rprospero@…> (01/25/18 12:05:09)

git-committer:

GitHub <noreply@…> (01/25/18 12:05:09)

Message:

Merge pull request #61 from rprospero/better_model_info_parser

Ticket 1035 - Simple Solution

Location:

sasmodels

Files:

• core.py (modified) (2 diffs)
• mixture.py (modified) (3 diffs)
• compare.py (modified) (12 diffs)
• details.py (modified) (3 diffs)
• generate.py (modified) (12 diffs)
• gengauss.py (added)
• guyou.py (moved) (moved from explore/guyou.py)
• jitter.py (moved) (moved from explore/jitter.py) (5 diffs)
• kernel_header.c (modified) (1 diff)
• kernel_iq.c (modified) (8 diffs)
• kernelpy.py (modified) (9 diffs)
• modelinfo.py (modified) (16 diffs)
• models/_spherepy.py (modified) (1 diff)
• models/barbell.c (modified) (3 diffs)
• models/bcc_paracrystal.c (modified) (3 diffs)
• models/capped_cylinder.c (modified) (4 diffs)
• models/core_shell_bicelle.c (modified) (2 diffs)
• models/core_shell_bicelle_elliptical.c (modified) (4 diffs)
• models/core_shell_bicelle_elliptical_belt_rough.c (modified) (6 diffs)
• models/core_shell_bicelle_elliptical_belt_rough.py (modified) (1 diff)
• models/core_shell_cylinder.c (modified) (3 diffs)
• models/core_shell_ellipsoid.c (modified) (3 diffs)
• models/core_shell_parallelepiped.c (modified) (3 diffs)
• models/core_shell_parallelepiped.py (modified) (7 diffs)
• models/cylinder.c (modified) (2 diffs)
• models/ellipsoid.c (modified) (2 diffs)
• models/elliptical_cylinder.c (modified) (3 diffs)
• models/elliptical_cylinder.py (modified) (1 diff)
• models/fcc_paracrystal.c (modified) (3 diffs)
• models/flexible_cylinder_elliptical.c (modified) (1 diff)
• models/hollow_cylinder.c (modified) (2 diffs)
• models/hollow_rectangular_prism.c (modified) (5 diffs)
• models/hollow_rectangular_prism.py (modified) (3 diffs)
• models/hollow_rectangular_prism_thin_walls.c (modified) (4 diffs)
• models/lib/gauss150.c (modified) (3 diffs)
• models/lib/gauss20.c (modified) (1 diff)
• models/lib/gauss76.c (modified) (2 diffs)
• models/line.py (modified) (2 diffs)
• models/parallelepiped.c (modified) (3 diffs)
• models/pringle.c (modified) (3 diffs)
• models/rectangular_prism.c (modified) (5 diffs)
• models/rectangular_prism.py (modified) (4 diffs)
• models/sc_paracrystal.c (modified) (3 diffs)
• models/stacked_disks.c (modified) (3 diffs)
• models/triaxial_ellipsoid.c (modified) (3 diffs)
• special.py (modified) (9 diffs)

sasmodels/core.py

@@ -148 +148 @@
     used with functions in generate to build the docs or extract model info.
     """
-    if '@' in model_string:
-        parts = model_string.split('@')
-        if len(parts) != 2:
-            raise ValueError("Use P@S to apply a structure factor S to model P")
-        P_info, Q_info = [load_model_info(part) for part in parts]
-        return product.make_product_info(P_info, Q_info)
-
-    product_parts = []
-    addition_parts = []
-
-    addition_parts_names = model_string.split('+')
-    if len(addition_parts_names) >= 2:
-        addition_parts = [load_model_info(part) for part in addition_parts_names]
-    elif len(addition_parts_names) == 1:
-        product_parts_names = model_string.split('*')
-        if len(product_parts_names) >= 2:
-            product_parts = [load_model_info(part) for part in product_parts_names]
-        elif len(product_parts_names) == 1:
-            if "custom." in product_parts_names[0]:
-                # Extract ModelName from "custom.ModelName"
-                pattern = "custom.([A-Za-z0-9_-]+)"
-                result = re.match(pattern, product_parts_names[0])
-                if result is None:
-                    raise ValueError("Model name in invalid format: " + product_parts_names[0])
-                model_name = result.group(1)
-                # Use ModelName to find the path to the custom model file
-                model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py")
-                if not os.path.isfile(model_path):
-                    raise ValueError("The model file {} doesn't exist".format(model_path))
-                kernel_module = custom.load_custom_kernel_module(model_path)
-                return modelinfo.make_model_info(kernel_module)
-            # Model is a core model
-            kernel_module = generate.load_kernel_module(product_parts_names[0])
-            return modelinfo.make_model_info(kernel_module)
-
-    model = None
-    if len(product_parts) > 1:
-        model = mixture.make_mixture_info(product_parts, operation='*')
-    if len(addition_parts) > 1:
-        if model is not None:
-            addition_parts.append(model)
-        model = mixture.make_mixture_info(addition_parts, operation='+')
-    return model
+    if "+" in model_string:
+        parts = [load_model_info(part)
+                 for part in model_string.split("+")]
+        return mixture.make_mixture_info(parts, operation='+')
+    elif "*" in model_string:
+        parts = [load_model_info(part)
+                 for part in model_string.split("*")]
+        return mixture.make_mixture_info(parts, operation='*')
+    elif "@" in model_string:
+        p_info, q_info = [load_model_info(part)
+                          for part in model_string.split("@")]
+        return product.make_product_info(p_info, q_info)
+    # We are now dealing with a pure model
+    elif "custom." in model_string:
+        pattern = "custom.([A-Za-z0-9_-]+)"
+        result = re.match(pattern, model_string)
+        if result is None:
+            raise ValueError("Model name in invalid format: " + model_string)
+        model_name = result.group(1)
+        # Use ModelName to find the path to the custom model file
+        model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py")
+        if not os.path.isfile(model_path):
+            raise ValueError("The model file {} doesn't exist".format(model_path))
+        kernel_module = custom.load_custom_kernel_module(model_path)
+        return modelinfo.make_model_info(kernel_module)
+    kernel_module = generate.load_kernel_module(model_string)
+    return modelinfo.make_model_info(kernel_module)
 
 
@@ -336 +320 @@
     print("\n".join(list_models(kind)))
 
+def test_load():
+    # type: () -> None
+    """Check that model load works"""
+    def test_models(fst, snd):
+        """Confirm that two models produce the same parameters"""
+        fst = load_model(fst)
+        snd = load_model(snd)
+        # Remove the upper case characters frin the parameters, since
+        # they lasndel the order of events and we specfically are
+        # changin that aspect
+        fst = [[x for x in p.name if x == x.lower()] for p in fst.info.parameters.kernel_parameters]
+        snd = [[x for x in p.name if x == x.lower()] for p in snd.info.parameters.kernel_parameters]
+        assert sorted(fst) == sorted(snd), "{} != {}".format(fst, snd)
+
+
+    test_models(
+        "cylinder+sphere",
+        "sphere+cylinder")
+    test_models(
+        "cylinder*sphere",
+        "sphere*cylinder")
+    test_models(
+        "cylinder@hardsphere*sphere",
+        "sphere*cylinder@hardsphere")
+    test_models(
+        "barbell+sphere*cylinder@hardsphere",
+        "sphere*cylinder@hardsphere+barbell")
+    test_models(
+        "barbell+cylinder@hardsphere*sphere",
+        "cylinder@hardsphere*sphere+barbell")
+    test_models(
+        "barbell+sphere*cylinder@hardsphere",
+        "barbell+cylinder@hardsphere*sphere")
+    test_models(
+        "sphere*cylinder@hardsphere+barbell",
+        "cylinder@hardsphere*sphere+barbell")
+    test_models(
+        "barbell+sphere*cylinder@hardsphere",
+        "cylinder@hardsphere*sphere+barbell")
+    test_models(
+        "barbell+cylinder@hardsphere*sphere",
+        "sphere*cylinder@hardsphere+barbell")
+
+    #Test the the model produces the parameters that we would expect
+    model = load_model("cylinder@hardsphere*sphere")
+    actual = [p.name for p in model.info.parameters.kernel_parameters]
+    target = ("sld sld_solvent radius length theta phi volfraction"
+              " A_sld A_sld_solvent A_radius").split()
+    assert target == actual, "%s != %s"%(target, actual)
+
 if __name__ == "__main__":
     list_models_main()

sasmodels/mixture.py

@@ -94 +94 @@
             # If model is a sum model, each constituent model gets its own scale parameter
             scale_prefix = prefix
-            if prefix == '' and part.operation == '*':
+            if prefix == '' and hasattr(part,"operation") and part.operation == '*':
                 # `part` is a composition product model. Find the prefixes of
                 # it's parameters to form a new prefix for the scale, eg:
@@ -233 +233 @@
         return self
 
-    def next(self):
+    def __next__(self):
         # type: () -> Tuple[List[Callable], CallDetails, np.ndarray]
         if self.part_num >= len(self.parts):
@@ -251 +251 @@
 
         return kernel, call_details, values
+
+    # CRUFT: py2 support
+    next = __next__
 
     def _part_details(self, info, par_index):

sasmodels/compare.py

@@ -42 +42 @@
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel, get_mesh
-from .generate import FLOAT_RE
+from .generate import FLOAT_RE, set_integration_size
 from .weights import plot_weights
 
@@ -92 +92 @@
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
     -neval=1 sets the number of evals for more accurate timing
+    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature
 
     === precision options ===
@@ -693 +694 @@
         data = empty_data2D(q, resolution=res)
         data.accuracy = opts['accuracy']
-        set_beam_stop(data, 0.0004)
+        set_beam_stop(data, qmin)
         index = ~data.mask
     else:
@@ -706 +707 @@
     return data, index
 
-def make_engine(model_info, data, dtype, cutoff):
+def make_engine(model_info, data, dtype, cutoff, ngauss=0):
     # type: (ModelInfo, Data, str, float) -> Calculator
     """
@@ -714 +715 @@
     than OpenCL.
     """
+    if ngauss:
+        set_integration_size(model_info, ngauss)
+
     if dtype is None or not dtype.endswith('!'):
         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
@@ -954 +958 @@
     'poly', 'mono', 'cutoff=',
     'magnetic', 'nonmagnetic',
-    'accuracy=',
+    'accuracy=', 'ngauss=',
     'neval=',  # for timing...
 
@@ -1089 +1093 @@
         'show_weights' : False,
         'sphere'    : 0,
+        'ngauss'    : '0',
     }
     for arg in flags:
@@ -1115 +1120 @@
         elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
         elif arg.startswith('-neval='):    opts['count'] = arg[7:]
+        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:]
         elif arg.startswith('-random='):
             opts['seed'] = int(arg[8:])
@@ -1169 +1175 @@
 
     comparison = any(PAR_SPLIT in v for v in values)
+
     if PAR_SPLIT in name:
         names = name.split(PAR_SPLIT, 2)
@@ -1181 +1188 @@
         return None
 
+    if PAR_SPLIT in opts['ngauss']:
+        opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        opts['ngauss'] = [int(opts['ngauss'])]*2
+
     if PAR_SPLIT in opts['engine']:
         opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
@@ -1199 +1212 @@
         opts['cutoff'] = [float(opts['cutoff'])]*2
 
-    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0])
+    base = make_engine(model_info[0], data, opts['engine'][0],
+                       opts['cutoff'][0], opts['ngauss'][0])
     if comparison:
-        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1])
+        comp = make_engine(model_info[1], data, opts['engine'][1],
+                           opts['cutoff'][1], opts['ngauss'][1])
     else:
         comp = None
@@ -1274 +1289 @@
         if model_info != model_info2:
             pars2 = randomize_pars(model_info2, pars2)
-            limit_dimensions(model_info, pars2, maxdim)
+            limit_dimensions(model_info2, pars2, maxdim)
             # Share values for parameters with the same name
             for k, v in pars.items():

sasmodels/details.py

@@ -258 +258 @@
     # type: (...) -> Sequence[np.ndarray]
     """
+    **Deprecated** Theta weights will be computed in the kernel wrapper if
+    they are needed.
+
     If there is a theta parameter, update the weights of that parameter so that
     the cosine weighting required for polar integration is preserved.
@@ -272 +275 @@
     Returns updated weights vectors
     """
-    # TODO: explain in a comment why scale and background are missing
     # Apparently the parameters.theta_offset similarly skips scale and
     # and background, so the indexing works out, but they are still shipped
@@ -279 +281 @@
         index = parameters.theta_offset
         theta = dispersity[index]
-        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,...
         theta_weight = abs(cos(radians(theta)))
         weights = tuple(theta_weight*w if k == index else w

sasmodels/generate.py

@@ -7 +7 @@
     particular dimensions averaged over all orientations.
 
-    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form
-    with particular dimensions for a single orientation.
-
-    *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the
-    polarized neutron spin states (up-up, up-down, down-up, down-down) for
-    a form with particular dimensions for a single orientation.
+    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc
+    for a rotationally symmetric form with particular dimensions.
+    qab, qc are determined from shape orientation and scattering angles.
+    This call is used if the shape has orientation parameters theta and phi.
+
+    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc
+    for a form with particular dimensions.  qa, qb, qc are determined from
+    shape orientation and scattering angles. This call is used if the shape
+    has orientation parameters theta, phi and psi.
+
+    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this
+    to create an arbitrary 2D theory function, needed for q-dependent
+    background functions and for models with non-uniform magnetism.
 
     *form_volume(p1, p2, ...)* returns the volume of the form with particular
@@ -31 +38 @@
 scale and background parameters for each model.
 
-*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99
-functions written for OpenCL.  All functions need prototype declarations
-even if the are defined before they are used.  OpenCL does not support
-*#include* preprocessor directives, so instead the list of includes needs
-to be given as part of the metadata in the kernel module definition.
-The included files should be listed using a path relative to the kernel
-module, or if using "lib/file.c" if it is one of the standard includes
-provided with the sasmodels source.  The includes need to be listed in
-order so that functions are defined before they are used.
+C code should be stylized C-99 functions written for OpenCL. All functions
+need prototype declarations even if the are defined before they are used.
+Although OpenCL supports *#include* preprocessor directives, the list of
+includes should be given as part of the metadata in the kernel module
+definition. The included files should be listed using a path relative to the
+kernel module, or if using "lib/file.c" if it is one of the standard includes
+provided with the sasmodels source. The includes need to be listed in order
+so that functions are defined before they are used.
 
 Floating point values should be declared as *double*.  For single precision
@@ -107 +113 @@
     present, the volume ratio is 1.
 
-    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the
-    C source code for the body of the volume, Iq, and Iqxy functions
+    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing
+    the C source code for the body of the volume, Iq, and Iqac functions
     respectively.  These can also be defined in the last source file.
 
-    *Iq* and *Iqxy* also be instead be python functions defining the
+    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the
     kernel.  If they are marked as *Iq.vectorized = True* then the
     kernel is passed the entire *q* vector at once, otherwise it is
@@ -168 +174 @@
 from zlib import crc32
 from inspect import currentframe, getframeinfo
+import logging
 
 import numpy as np  # type: ignore
@@ -181 +188 @@
     pass
 # pylint: enable=unused-import
+
+logger = logging.getLogger(__name__)
 
 # jitter projection to use in the kernel code.  See explore/jitter.py
@@ -270 +279 @@
 """
 
+
+def set_integration_size(info, n):
+    # type: (ModelInfo, int) -> None
+    """
+    Update the model definition, replacing the gaussian integration with
+    a gaussian integration of a different size.
+
+    Note: this really ought to be a method in modelinfo, but that leads to
+    import loops.
+    """
+    if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)):
+        import os.path
+        from .gengauss import gengauss
+        path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n)
+        if not os.path.exists(path):
+            gengauss(n, path)
+        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
+                        else lib for lib in info.source]
 
 def format_units(units):
@@ -608 +635 @@
 
 """
-def _gen_fn(name, pars, body, filename, line):
-    # type: (str, List[Parameter], str, str, int) -> str
+def _gen_fn(model_info, name, pars):
+    # type: (ModelInfo, str, List[Parameter]) -> str
     """
     Generate a function given pars and body.
@@ -621 +648 @@
     """
     par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
+    body = getattr(model_info, name)
+    filename = model_info.filename
+    # Note: if symbol is defined strangely in the module then default it to 1
+    lineno = model_info.lineno.get(name, 1)
     return _FN_TEMPLATE % {
         'name': name, 'pars': par_decl, 'body': body,
-        'filename': filename.replace('\\', '\\\\'), 'line': line,
+        'filename': filename.replace('\\', '\\\\'), 'line': lineno,
     }
 
@@ -638 +669 @@
 
 # type in IQXY pattern could be single, float, double, long double, ...
-_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)",
+_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]",
                            flags=re.MULTILINE)
-def _have_Iqxy(sources):
+def find_xy_mode(source):
     # type: (List[str]) -> bool
     """
-    Return true if any file defines Iqxy.
+    Return the xy mode as qa, qac, qabc or qxy.
 
     Note this is not a C parser, and so can be easily confused by
     non-standard syntax.  Also, it will incorrectly identify the following
-    as having Iqxy::
+    as having 2D models::
 
         /*
-        double Iqxy(qx, qy, ...) { ... fill this in later ... }
+        double Iqac(qab, qc, ...) { ... fill this in later ... }
         */
 
-    If you want to comment out an Iqxy function, use // on the front of the
-    line instead.
-    """
-    for _path, code in sources:
-        if _IQXY_PATTERN.search(code):
-            return True
-    return False
-
-
-def _add_source(source, code, path):
+    If you want to comment out the function, use // on the front of the
+    line::
+
+        /*
+        // double Iqac(qab, qc, ...) { ... fill this in later ... }
+        */
+
+    """
+    for code in source:
+        m = _IQXY_PATTERN.search(code)
+        if m is not None:
+            return m.group('mode')
+    return 'qa'
+
+
+def _add_source(source, code, path, lineno=1):
     """
     Add a file to the list of source code chunks, tagged with path and line.
     """
     path = path.replace('\\', '\\\\')
-    source.append('#line 1 "%s"' % path)
+    source.append('#line %d "%s"' % (lineno, path))
     source.append(code)
 
@@ -698 +735 @@
     user_code = [(f, open(f).read()) for f in model_sources(model_info)]
 
-    # What kind of 2D model do we need?
-    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str)
-               else 'qac' if not partable.is_asymmetric
-               else 'qabc')
-
     # Build initial sources
     source = []
@@ -709 +741 @@
         _add_source(source, code, path)
 
+    if model_info.c_code:
+        _add_source(source, model_info.c_code, model_info.filename,
+                    lineno=model_info.lineno.get('c_code', 1))
+
     # Make parameters for q, qx, qy so that we can use them in declarations
-    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
+    q, qx, qy, qab, qa, qb, qc \
+        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
     # Generate form_volume function, etc. from body only
     if isinstance(model_info.form_volume, str):
         pars = partable.form_volume_parameters
-        source.append(_gen_fn('form_volume', pars, model_info.form_volume,
-                              model_info.filename, model_info._form_volume_line))
+        source.append(_gen_fn(model_info, 'form_volume', pars))
     if isinstance(model_info.Iq, str):
         pars = [q] + partable.iq_parameters
-        source.append(_gen_fn('Iq', pars, model_info.Iq,
-                              model_info.filename, model_info._Iq_line))
+        source.append(_gen_fn(model_info, 'Iq', pars))
     if isinstance(model_info.Iqxy, str):
-        pars = [qx, qy] + partable.iqxy_parameters
-        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy,
-                              model_info.filename, model_info._Iqxy_line))
+        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters
+        source.append(_gen_fn(model_info, 'Iqxy', pars))
+    if isinstance(model_info.Iqac, str):
+        pars = [qab, qc] + partable.iq_parameters
+        source.append(_gen_fn(model_info, 'Iqac', pars))
+    if isinstance(model_info.Iqabc, str):
+        pars = [qa, qb, qc] + partable.iq_parameters
+        source.append(_gen_fn(model_info, 'Iqabc', pars))
+
+    # What kind of 2D model do we need?  Is it consistent with the parameters?
+    xy_mode = find_xy_mode(source)
+    if xy_mode == 'qabc' and not partable.is_asymmetric:
+        raise ValueError("asymmetric oriented models need to define Iqabc")
+    elif xy_mode == 'qac' and partable.is_asymmetric:
+        raise ValueError("symmetric oriented models need to define Iqac")
+    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'):
+        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode)
+    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'):
+        if xy_mode == 'qxy':
+            logger.warn("oriented shapes should define Iqac or Iqabc")
+        else:
+            raise ValueError("Expected function Iqac or Iqabc for oriented shape")
 
     # Define the parameter table
@@ -749 +803 @@
     if xy_mode == 'qabc':
         pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
-        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars
+        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_ABC"
     elif xy_mode == 'qac':
         pars = ",".join(["_qa", "_qc"] + model_refs)
-        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars
+        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_AC"
-    else:  # xy_mode == 'qa'
+    elif xy_mode == 'qa':
         pars = ",".join(["_qa"] + model_refs)
         call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_A"
+    elif xy_mode == 'qxy':
+        orientation_refs = _call_pars("_v.", partable.orientation_parameters)
+        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs)
+        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars
+        clear_iqxy = "#undef CALL_IQ_XY"
+        if partable.orientation_parameters:
+            call_iqxy += "\n#define HAVE_THETA"
+            clear_iqxy += "\n#undef HAVE_THETA"
+        if partable.is_asymmetric:
+            call_iqxy += "\n#define HAVE_PSI"
+            clear_iqxy += "\n#undef HAVE_PSI"
+
 
     magpars = [k-2 for k, p in enumerate(partable.call_parameters)

sasmodels/jitter.py

@@ -1 +1 @@
 #!/usr/bin/env python
 """
-Application to explore the difference between sasview 3.x orientation
-dispersity and possible replacement algorithms.
+Jitter Explorer
+===============
+
+Application to explore orientation angle and angular dispersity.
 """
 from __future__ import division, print_function
@@ -124 +126 @@
 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
     """Draw a parallelepiped."""
-    a,b,c = size
+    a, b, c = size
     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
@@ -142 +144 @@
     ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
 
+    # Draw pink face on box.
+    # Since I can't control face color, instead draw a thin box situated just
+    # in front of the "a+" face.
+    if 1:
+        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
+        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
+        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
+        x, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z)
+        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha)
+
     draw_labels(ax, view, jitter, [
          ('c+', [ 0, 0, c], [ 1, 0, 0]),
@@ -165 +177 @@
     # constants in kernel_iq.c
     'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
+    'azimuthal_equal_area',
 ]
 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
@@ -607 +620 @@
     Show an interactive orientation and jitter demo.
 
-    *model_name* is one of the models available in :func:`select_model`.
+    *model_name* is one of: sphere, cylinder, ellipsoid, parallelepiped,
+    triaxial_ellipsoid, or bcc_paracrystal.
+
+    *size* gives the dimensions (a, b, c) of the shape.
+
+    *dist* is the type of dispersition: gaussian, rectangle, or uniform.
+
+    *mesh* is the number of points in the dispersion mesh.
+
+    *projection* is the map projection to use for the mesh: equirectangular,
+    sinusoidal, guyou, azimuthal_equidistance, or azimuthal_equal_area.
     """
     # projection number according to 1-order position in list, but

sasmodels/kernel_header.c

@@ -150 +150 @@
 inline double cube(double x) { return x*x*x; }
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
+
+// CRUFT: support old style models with orientation received qx, qy and angles
+
+// To rotate from the canonical position to theta, phi, psi, first rotate by
+// psi about the major axis, oriented along z, which is a rotation in the
+// detector plane xy. Next rotate by theta about the y axis, aligning the major
+// axis in the xz plane. Finally, rotate by phi in the detector plane xy.
+// To compute the scattering, undo these rotations in reverse order:
+//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi
+// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit
+// vector in the q direction.
+// To change between counterclockwise and clockwise rotation, change the
+// sign of phi and psi.
+
+#if 1
+//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017
+#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
+    SINCOS(phi*M_PI_180, sn, cn); \
+    q = sqrt(qx*qx + qy*qy); \
+    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \
+    sn = sqrt(1 - cn*cn); \
+    } while (0)
+#else
+// SasView 3.x definition of orientation
+#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
+    SINCOS(theta*M_PI_180, sn, cn); \
+    q = sqrt(qx*qx + qy*qy);\
+    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \
+    sn = sqrt(1 - cn*cn); \
+    } while (0)
+#endif
+
+#if 1
+#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \
+    q = sqrt(qx*qx + qy*qy); \
+    const double qxhat = qx/q; \
+    const double qyhat = qy/q; \
+    double sin_theta, cos_theta; \
+    double sin_phi, cos_phi; \
+    double sin_psi, cos_psi; \
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
+    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \
+         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \
+    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \
+         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \
+    zhat = qxhat*(-sin_theta*cos_phi) \
+         + qyhat*(-sin_theta*sin_phi); \
+    } while (0)
+#else
+// SasView 3.x definition of orientation
+#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \
+    q = sqrt(qx*qx + qy*qy); \
+    const double qxhat = qx/q; \
+    const double qyhat = qy/q; \
+    double sin_theta, cos_theta; \
+    double sin_phi, cos_phi; \
+    double sin_psi, cos_psi; \
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
+    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \
+    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \
+    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \
+    } while (0)
+#endif

sasmodels/kernel_iq.c

@@ -31 +31 @@
 //  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
 //  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
+//  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models
 //  INVALID(table) : test if the current point is feesible to calculate.  This
 //      will be defined in the kernel definition file.
@@ -173 +174 @@
 static double
 qac_apply(
-    QACRotation rotation,
+    QACRotation *rotation,
     double qx, double qy,
     double *qa_out, double *qc_out)
 {
-    const double dqc = rotation.R31*qx + rotation.R32*qy;
+    const double dqc = rotation->R31*qx + rotation->R32*qy;
     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
@@ -246 +247 @@
 static double
 qabc_apply(
-    QABCRotation rotation,
+    QABCRotation *rotation,
     double qx, double qy,
     double *qa_out, double *qb_out, double *qc_out)
 {
-    *qa_out = rotation.R11*qx + rotation.R12*qy;
-    *qb_out = rotation.R21*qx + rotation.R22*qy;
-    *qc_out = rotation.R31*qx + rotation.R32*qy;
+    *qa_out = rotation->R11*qx + rotation->R12*qy;
+    *qb_out = rotation->R21*qx + rotation->R22*qy;
+    *qc_out = rotation->R31*qx + rotation->R32*qy;
 }
 
@@ -453 +454 @@
   // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
   #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
-  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
+  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)
   #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
 
@@ -467 +468 @@
   local_values.table.psi = 0.;
   #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
-  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
+  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
-#endif
-
-// Doing jitter projection code outside the previous if block so that we don't
-// need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This
-// will become more important if we implement more projections, or more
-// complicated projections.
-#if defined(CALL_IQ) || defined(CALL_IQ_A)
+#elif defined(CALL_IQ_XY)
+  // direct call to qx,qy calculator
+  double qx, qy;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table)
+#endif
+
+// Define APPLY_PROJECTION depending on model symmetries. We do this outside
+// the previous if block so that we don't need to repeat the identical
+// logic in the IQ_AC and IQ_ABC branches.  This will become more important
+// if we implement more projections, or more complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
-#else // !spherosymmetric projection
+#elif defined(CALL_IQ_XY) // pass orientation to the model
+  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc
+  // Need to plug the values for the orientation angles back into parameter
+  // table in case they were overridden by the orientation offset.  This
+  // means that orientation dispersity will not work for these models, but
+  // it was broken anyway, so no matter.  Still want to provide Iqxy in case
+  // the user model wants full control of orientation/magnetism.
+  #if defined(HAVE_PSI)
+    const double theta = values[details->theta_par+2];
+    const double phi = values[details->theta_par+3];
+    const double psi = values[details->theta_par+4];
+    double weight;
+    #define APPLY_PROJECTION() do { \
+      local_values.table.theta = theta; \
+      local_values.table.phi = phi; \
+      local_values.table.psi = psi; \
+      weight=weight0; \
+    } while (0)
+  #elif defined(HAVE_THETA)
+    const double theta = values[details->theta_par+2];
+    const double phi = values[details->theta_par+3];
+    double weight;
+    #define APPLY_PROJECTION() do { \
+      local_values.table.theta = theta; \
+      local_values.table.phi = phi; \
+      weight=weight0; \
+    } while (0)
+  #else
+    #define APPLY_PROJECTION() const double weight=weight0
+  #endif
+#else // apply jitter and view before calling the model
   // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
   const double theta = values[details->theta_par+2];
@@ -488 +526 @@
   // we go through the mesh.
   double dtheta, dphi, weight;
-  #if PROJECTION == 1
+  #if PROJECTION == 1 // equirectangular
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -494 +532 @@
       weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
     } while (0)
-  #elif PROJECTION == 2
+  #elif PROJECTION == 2 // sinusoidal
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
@@ -504 +542 @@
     } while (0)
   #endif
-#endif // !spherosymmetric projection
+#endif // done defining APPLY_PROJECTION
 
 // ** define looping macros **

sasmodels/kernelpy.py

@@ -26 +26 @@
 # pylint: enable=unused-import
 
+logger = logging.getLogger(__name__)
+
 class PyModel(KernelModel):
     """
@@ -31 +33 @@
     """
     def __init__(self, model_info):
-        # Make sure Iq and Iqxy are available and vectorized
+        # Make sure Iq is available and vectorized
         _create_default_functions(model_info)
         self.info = model_info
         self.dtype = np.dtype('d')
-        logging.info("load python model " + self.info.name)
+        logger.info("load python model " + self.info.name)
 
     def make_kernel(self, q_vectors):
         q_input = PyInput(q_vectors, dtype=F64)
-        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq
-        return PyKernel(kernel, self.info, q_input)
+        return PyKernel(self.info, q_input)
 
     def release(self):
@@ -89 +90 @@
     Callable SAS kernel.
 
-    *kernel* is the DllKernel object to call.
+    *kernel* is the kernel object to call.
 
     *model_info* is the module information
@@ -104 +105 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, model_info, q_input):
+    def __init__(self, model_info, q_input):
         # type: (callable, ModelInfo, List[np.ndarray]) -> None
         self.dtype = np.dtype('d')
@@ -110 +111 @@
         self.q_input = q_input
         self.res = np.empty(q_input.nq, q_input.dtype)
-        self.kernel = kernel
         self.dim = '2d' if q_input.is_2d else '1d'
 
@@ -159 +159 @@
         # Generate a closure which calls the form_volume if it exists.
         form_volume = model_info.form_volume
-        self._volume = ((lambda: form_volume(*volume_args)) if form_volume
-                        else (lambda: 1.0))
+        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else
+                        (lambda: 1.0))
 
     def __call__(self, call_details, values, cutoff, magnetic):
@@ -261 +261 @@
     any functions that are not already marked as vectorized.
     """
+    # Note: must call create_vector_Iq before create_vector_Iqxy
     _create_vector_Iq(model_info)
-    _create_vector_Iqxy(model_info)  # call create_vector_Iq() first
+    _create_vector_Iqxy(model_info)
 
 
@@ -280 +281 @@
         model_info.Iq = vector_Iq
 
+
 def _create_vector_Iqxy(model_info):
     """
     Define Iqxy as a vector function if it exists, or default it from Iq().
     """
-    Iq, Iqxy = model_info.Iq, model_info.Iqxy
+    Iqxy = getattr(model_info, 'Iqxy', None)
     if callable(Iqxy):
         if not getattr(Iqxy, 'vectorized', False):
@@ -295 +297 @@
             vector_Iqxy.vectorized = True
             model_info.Iqxy = vector_Iqxy
-    elif callable(Iq):
+    else:
         #print("defaulting Iqxy")
         # Iq is vectorized because create_vector_Iq was already called.
+        Iq = model_info.Iq
         def default_Iqxy(qx, qy, *args):
             """

sasmodels/modelinfo.py

@@ -37 +37 @@
 
 # assumptions about common parameters exist throughout the code, such as:
-# (1) kernel functions Iq, Iqxy, form_volume, ... don't see them
+# (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them
 # (2) kernel drivers assume scale is par[0] and background is par[1]
 # (3) mixture models drop the background on components and replace the scale
@@ -256 +256 @@
 
     *type* indicates how the parameter will be used.  "volume" parameters
-    will be used in all functions.  "orientation" parameters will be used
-    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
-    *Imagnetic* only.  If *type* is the empty string, the parameter will
+    will be used in all functions.  "orientation" parameters are not passed,
+    but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to
+    *Iqxy* and *Imagnetic*.  If *type* is the empty string, the parameter will
     be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters
     can automatically be promoted to magnetic parameters, each of which
@@ -386 +386 @@
       with vector parameter p sent as p[].
 
-    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
-      function, with vector parameter p sent as p[].
-
     * *form_volume_parameters* is the list of parameters to the form_volume(...)
       function, with vector parameter p sent as p[].
@@ -443 +440 @@
         self.iq_parameters = [p for p in self.kernel_parameters
                               if p.type not in ('orientation', 'magnetic')]
-        # note: orientation no longer sent to Iqxy, so its the same as
-        #self.iqxy_parameters = [p for p in self.kernel_parameters
-        #                        if p.type != 'magnetic']
+        self.orientation_parameters = [p for p in self.kernel_parameters
+                                       if p.type == 'orientation']
         self.form_volume_parameters = [p for p in self.kernel_parameters
                                        if p.type == 'volume']
@@ -490 +486 @@
                 if p.type != 'orientation':
                     raise TypeError("psi must be an orientation parameter")
+            elif p.type == 'orientation':
+                raise TypeError("only theta, phi and psi can be orientation parameters")
         if theta >= 0 and phi >= 0:
+            last_par = len(self.kernel_parameters) - 1
             if phi != theta+1:
                 raise TypeError("phi must follow theta")
             if psi >= 0 and psi != phi+1:
                 raise TypeError("psi must follow phi")
+            #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par):
+            #    raise TypeError("orientation parameters must appear at the "
+            #                    "end of the parameter table")
         elif theta >= 0 or phi >= 0 or psi >= 0:
             raise TypeError("oriented shapes must have both theta and phi and maybe psi")
@@ -715 +717 @@
 
 
+#: Set of variables defined in the model that might contain C code
+C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code']
+
 def _find_source_lines(model_info, kernel_module):
     # type: (ModelInfo, ModuleType) -> None
@@ -720 +725 @@
     Identify the location of the C source inside the model definition file.
 
-    This code runs through the source of the kernel module looking for
-    lines that start with 'Iq', 'Iqxy' or 'form_volume'.  Clearly there are
-    all sorts of reasons why this might not work (e.g., code commented out
-    in a triple-quoted line block, code built using string concatenation,
-    or code defined in the branch of an 'if' block), but it should work
-    properly in the 95% case, and getting the incorrect line number will
-    be harmless.
-    """
-    # Check if we need line numbers at all
-    if callable(model_info.Iq):
-        return None
-
-    if (model_info.Iq is None
-            and model_info.Iqxy is None
-            and model_info.Imagnetic is None
-            and model_info.form_volume is None):
+    This code runs through the source of the kernel module looking for lines
+    that contain C code (because it is a c function definition).  Clearly
+    there are all sorts of reasons why this might not work (e.g., code
+    commented out in a triple-quoted line block, code built using string
+    concatenation, code defined in the branch of an 'if' block, code imported
+    from another file), but it should work properly in the 95% case, and for
+    the remainder, getting the incorrect line number will merely be
+    inconvenient.
+    """
+    # Only need line numbers if we are creating a C module and the C symbols
+    # are defined.
+    if (callable(model_info.Iq)
+            or not any(hasattr(model_info, s) for s in C_SYMBOLS)):
         return
 
-    # find the defintion lines for the different code blocks
+    # load the module source if we can
     try:
         source = inspect.getsource(kernel_module)
     except IOError:
         return
-    for k, v in enumerate(source.split('\n')):
-        if v.startswith('Imagnetic'):
-            model_info._Imagnetic_line = k+1
-        elif v.startswith('Iqxy'):
-            model_info._Iqxy_line = k+1
-        elif v.startswith('Iq'):
-            model_info._Iq_line = k+1
-        elif v.startswith('form_volume'):
-            model_info._form_volume_line = k+1
-
+
+    # look for symbol at the start of the line
+    for lineno, line in enumerate(source.split('\n')):
+        for name in C_SYMBOLS:
+            if line.startswith(name):
+                # Add 1 since some compilers complain about "#line 0"
+                model_info.lineno[name] = lineno + 1
+                break
 
 def make_model_info(kernel_module):
@@ -761 +761 @@
     Fill in default values for parts of the module that are not provided.
 
-    Note: vectorized Iq and Iqxy functions will be created for python
+    Note: vectorized Iq and Iqac/Iqabc functions will be created for python
     models when the model is first called, not when the model is loaded.
     """
@@ -790 +790 @@
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
     info.source = getattr(kernel_module, 'source', [])
+    info.c_code = getattr(kernel_module, 'c_code', None)
     # TODO: check the structure of the tests
     info.tests = getattr(kernel_module, 'tests', [])
@@ -797 +798 @@
     info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
     info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
+    info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore
+    info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore
     info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore
     info.profile = getattr(kernel_module, 'profile', None) # type: ignore
@@ -811 +814 @@
     info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore
 
+    if callable(info.Iq) and parameters.has_2d:
+        raise ValueError("oriented python models not supported")
+
+    info.lineno = {}
     _find_source_lines(info, kernel_module)
 
@@ -824 +831 @@
 
     The structure should be mostly static, other than the delayed definition
-    of *Iq* and *Iqxy* if they need to be defined.
+    of *Iq*, *Iqac* and *Iqabc* if they need to be defined.
     """
     #: Full path to the file defining the kernel, if any.
@@ -906 +913 @@
     structure_factor = None # type: bool
     #: List of C source files used to define the model.  The source files
-    #: should define the *Iq* function, and possibly *Iqxy*, though a default
-    #: *Iqxy = Iq(sqrt(qx**2+qy**2)* will be created if no *Iqxy* is provided.
-    #: Files containing the most basic functions must appear first in the list,
-    #: followed by the files that use those functions.  Form factors are
-    #: indicated by providing a :attr:`ER` function.
+    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
+    #: model defines orientation parameters. Files containing the most basic
+    #: functions must appear first in the list, followed by the files that
+    #: use those functions.  Form factors are indicated by providing
+    #: an :attr:`ER` function.
     source = None           # type: List[str]
     #: The set of tests that must pass.  The format of the tests is described
@@ -935 +942 @@
     #: See :attr:`ER` for details on the parameters.
     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]]
+    #: Arbitrary C code containing supporting functions, etc., to be inserted
+    #: after everything in source.  This can include Iq and Iqxy functions with
+    #: the full function signature, including all parameters.
+    c_code = None
     #: Returns the form volume for python-based models.  Form volume is needed
     #: for volume normalization in the polydispersity integral.  If no
@@ -955 +966 @@
     #: include the decimal point. See :mod:`generate` for more details.
     Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
-    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
-    Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
+    Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
+    Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
     Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
@@ -972 +985 @@
     #: Returns a random parameter set for the model
     random = None           # type: Optional[Callable[[], Dict[str, float]]]
-
-    # line numbers within the python file for bits of C source, if defined
-    # NB: some compilers fail with a "#line 0" directive, so default to 1.
-    _Imagnetic_line = 1
-    _Iqxy_line = 1
-    _Iq_line = 1
-    _form_volume_line = 1
-
+    #: Line numbers for symbols defining C code
+    lineno = None           # type: Dict[str, int]
 
     def __init__(self):

sasmodels/models/_spherepy.py

@@ -88 +88 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
-def Iqxy(qx, qy, sld, sld_solvent, radius):
-    return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius)
-Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
-
 def sesans(z, sld, sld_solvent, radius):
     """

sasmodels/models/barbell.c

@@ -23 +23 @@
     const double qab_r = radius_bell*qab; // Q*R*sin(theta)
     double total = 0.0;
-    for (int i = 0; i < 76; i++){
-        const double t = Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++){
+        const double t = GAUSS_Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += Gauss76Wt[i] * Fq;
+        total += GAUSS_W[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -73 +73 @@
     const double zb = M_PI_4;
     double total = 0.0;
-    for (int i = 0; i < 76; i++){
-        const double alpha= Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++){
+        const double alpha= GAUSS_Z[i]*zm + zb;
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
         const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length);
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        total += GAUSS_W[i] * Aq * Aq * sin_alpha;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -90 +90 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld, double solvent_sld,
     double radius_bell, double radius, double length)

sasmodels/models/bcc_paracrystal.c

@@ -81 +81 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -95 +95 @@
             const double qb = qab*sin_phi;
             const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;
@@ -107 +107 @@
 
 
-static double Iqxy(double qa, double qb, double qc,
+static double Iqabc(double qa, double qb, double qc,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld)

sasmodels/models/capped_cylinder.c

@@ -30 +30 @@
     const double qab_r = radius_cap*qab; // Q*R*sin(theta)
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double t = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N; i++) {
+        const double t = GAUSS_Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += Gauss76Wt[i] * Fq;
+        total += GAUSS_W[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -95 +95 @@
     const double zb = M_PI_4;
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double theta = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N ;i++) {
+        const double theta = GAUSS_Z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
@@ -103 +103 @@
         const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
         // scale by sin_theta for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_theta;
+        total += GAUSS_W[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -115 +115 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld, double solvent_sld, double radius,
     double radius_cap, double length)

sasmodels/models/core_shell_bicelle.c

@@ -52 +52 @@
 
     double total = 0.0;
-    for(int i=0;i<N_POINTS_76;i++) {
-        double theta = (Gauss76Z[i] + 1.0)*uplim;
+    for(int i=0;i<GAUSS_N;i++) {
+        double theta = (GAUSS_Z[i] + 1.0)*uplim;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
                                    halflength, sld_core, sld_face, sld_rim, sld_solvent);
-        total += Gauss76Wt[i]*fq*fq*sin_theta;
+        total += GAUSS_W[i]*fq*fq*sin_theta;
     }
 
@@ -67 +67 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double radius,
     double thick_rim,

sasmodels/models/core_shell_bicelle_elliptical.c

@@ -37 +37 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         //const double va = 0.0;
         //const double vb = 1.0;
-        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0;
+        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0;
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         const double qab = q*sin_theta;
@@ -49 +49 @@
         const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
         double inner_sum=0.0;
-        for(int j=0;j<76;j++) {
+        for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
             // inner integral limits
             //const double vaj=0.0;
             //const double vbj=M_PI;
-            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2;
+            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(phi));
             const double be1 = sas_2J1x_x(rr*qab);
@@ -61 +61 @@
             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
 
-            inner_sum += Gauss76Wt[j] * fq * fq;
+            inner_sum += GAUSS_W[j] * fq * fq;
         }
         //now calculate outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
 
@@ -71 +71 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double r_minor,
     double x_core,

sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

@@ -7 +7 @@
         double length)
 {
-    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
+    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +
                  square(r_minor)*x_core*2.0*thick_face  );
 }
@@ -47 +47 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         // since we generate these lots of times, why not store them somewhere?
-        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
+        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0;
         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
         double inner_sum=0;
@@ -58 +58 @@
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
-        for(int j=0;j<76;j++) {
+        for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
-            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
+            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(beta));
             double besarg1 = q*rr*sin_alpha;
@@ -67 +67 @@
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
-            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
+            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 +
                                               dr2*si1*be2 +
                                               dr3*si2*be1);
         }
         //now calculate outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
 
@@ -79 +79 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
           double r_minor,
           double x_core,
@@ -114 +114 @@
     return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma));
 }
-

sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

@@ -149 +149 @@
     ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"],
     ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"],
+    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"],
     ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"],
     ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"],
     ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"],
-    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"]
     ]
 

sasmodels/models/core_shell_cylinder.c

@@ -30 +30 @@
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
+    for (int i=0; i<GAUSS_N ;i++) {
         // translate a point in [-1,1] to a point in [0, pi/2]
-        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0;
         double sin_theta, cos_theta;
-        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4;
+        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4;
         SINCOS(theta, sin_theta,  cos_theta);
         const double qab = q*sin_theta;
@@ -40 +40 @@
         const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
             + _cyl(shell_vd, shell_r*qab, shell_h*qc);
-        total += Gauss76Wt[i] * fq * fq * sin_theta;
+        total += GAUSS_W[i] * fq * fq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -48 +48 @@
 
 
-double Iqxy(double qab, double qc,
+double Iqac(double qab, double qc,
     double core_sld,
     double shell_sld,

sasmodels/models/core_shell_ellipsoid.c

@@ -59 +59 @@
     const double b = 0.5;
     double total = 0.0;     //initialize intergral
-    for(int i=0;i<76;i++) {
-        const double cos_theta = Gauss76Z[i]*m + b;
+    for(int i=0;i<GAUSS_N;i++) {
+        const double cos_theta = GAUSS_Z[i]*m + b;
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta,
@@ -66 +66 @@
             equat_shell, polar_shell,
             sld_core_shell, sld_shell_solvent);
-        total += Gauss76Wt[i] * fq * fq;
+        total += GAUSS_W[i] * fq * fq;
     }
     total *= m;
@@ -75 +75 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double radius_equat_core,
     double x_core,

sasmodels/models/core_shell_parallelepiped.c

@@ +1 @@
+// Set OVERLAPPING to 1 in order to fill in the edges of the box, with
+// c endcaps and b overlapping a.  With the proper choice of parameters,
+// (setting rim slds to sld, core sld to solvent, rim thickness to thickness
+// and subtracting 2*thickness from length, this should match the hollow
+// rectangular prism.)  Set it to 0 for the documented behaviour.
+#define OVERLAPPING 0
 static double
 form_volume(double length_a, double length_b, double length_c,
     double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
-    //return length_a * length_b * length_c;
-    return length_a * length_b * length_c +
-           2.0 * thick_rim_a * length_b * length_c +
-           2.0 * thick_rim_b * length_a * length_c +
-           2.0 * thick_rim_c * length_a * length_b;
+    return
+#if OVERLAPPING
+        // Hollow rectangular prism only includes the volume of the shell
+        // so uncomment the next line when comparing.  Solid rectangular
+        // prism, or parallelepiped want filled cores, so comment when
+        // comparing.
+        //-length_a * length_b * length_c +
+        (length_a + 2.0*thick_rim_a) *
+        (length_b + 2.0*thick_rim_b) *
+        (length_c + 2.0*thick_rim_c);
+#else
+        length_a * length_b * length_c +
+        2.0 * thick_rim_a * length_b * length_c +
+        2.0 * length_a * thick_rim_b * length_c +
+        2.0 * length_a * length_b * thick_rim_c;
+#endif
 }
 
@@ -24 +41 @@
     double thick_rim_c)
 {
-    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
+    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
+    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
+    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
 
-    const double mu = 0.5 * q * length_b;
+    const double half_q = 0.5*q;
 
-    //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
-    //double vol = length_a * length_b * length_c +
-    //       2.0 * thick_rim_a * length_b * length_c +
-    //       2.0 * thick_rim_b * length_a * length_c +
-    //       2.0 * thick_rim_c * length_a * length_b;
+    const double tA = length_a + 2.0*thick_rim_a;
+    const double tB = length_b + 2.0*thick_rim_b;
+    const double tC = length_c + 2.0*thick_rim_c;
 
-    // Scale sides by B
-    const double a_scaled = length_a / length_b;
-    const double c_scaled = length_c / length_b;
-
-    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
-    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
-    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
-
-    double Vin = length_a * length_b * length_c;
-    //double Vot = (length_a * length_b * length_c +
-    //            2.0 * thick_rim_a * length_b * length_c +
-    //            2.0 * length_a * thick_rim_b * length_c +
-    //            2.0 * length_a * length_b * thick_rim_c);
-    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
-    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
-
-    // Scale factors (note that drC is not used later)
-    const double drho0 = (core_sld-solvent_sld);
-    const double drhoA = (arim_sld-solvent_sld);
-    const double drhoB = (brim_sld-solvent_sld);
-    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
-
-
-    // Precompute scale factors for combining cross terms from the shape
-    const double scale23 = drhoA*V1;
-    const double scale14 = drhoB*V2;
-    const double scale24 = drhoC*V3;
-    const double scale11 = drho0*Vin;
-    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
+    // Scale factors
+    const double dr0 = (core_sld-solvent_sld);
+    const double drA = (arim_sld-solvent_sld);
+    const double drB = (brim_sld-solvent_sld);
+    const double drC = (crim_sld-solvent_sld);
 
     // outer integral (with gauss points), integration limits = 0, 1
-    double outer_total = 0; //initialize integral
+    double outer_sum = 0; //initialize integral
+    for( int i=0; i<GAUSS_N; i++) {
+        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
+        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
 
-    for( int i=0; i<76; i++) {
-        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
-        double mu_proj = mu * sqrt(1.0-sigma*sigma);
+        // inner integral (with gauss points), integration limits = 0, pi/2
+        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
+        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
+        double inner_sum = 0.0;
+        for(int j=0; j<GAUSS_N; j++) {
+            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 );
+            double sin_beta, cos_beta;
+            SINCOS(M_PI_2*beta, sin_beta, cos_beta);
+            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
+            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
+            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
+            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
 
-        // inner integral (with gauss points), integration limits = 0, 1
-        double inner_total = 0.0;
-        double inner_total_crim = 0.0;
-        for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
-            double sin_uu, cos_uu;
-            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
-            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
-            const double si2 = sas_sinx_x(mu_proj * cos_uu );
-            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
-            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
+#if OVERLAPPING
+            const double f = dr0*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siAt*(siBt-siB)*siC
+                + drC*siAt*siBt*(siCt-siC);
+#else
+            const double f = dr0*siA*siB*siC
+                + drA*(siAt-siA)*siB*siC
+                + drB*siA*(siBt-siB)*siC
+                + drC*siA*siB*(siCt-siC);
+#endif
 
-            // Expression in libCylinder.c (neither drC nor Vot are used)
-            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
-            const double form_crim = scale11*si1*si2;
-
-            //  correct FF : sum of square of phase factors
-            inner_total += Gauss76Wt[j] * form * form;
-            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
+            inner_sum += GAUSS_W[j] * f * f;
         }
-        inner_total *= 0.5;
-        inner_total_crim *= 0.5;
+        inner_sum *= 0.5;
         // now sum up the outer integral
-        const double si = sas_sinx_x(mu * c_scaled * sigma);
-        const double si_crim = sas_sinx_x(mu * tc * sigma);
-        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
-    outer_total *= 0.5;
+    outer_sum *= 0.5;
 
     //convert from [1e-12 A-1] to [cm-1]
-    return 1.0e-4 * outer_total;
+    return 1.0e-4 * outer_sum;
 }
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double core_sld,
     double arim_sld,
@@ -128 +121 @@
     const double drC = crim_sld-solvent_sld;
 
-    double Vin = length_a * length_b * length_c;
-    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc)
-    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = 2.0 * length_a * length_b * thick_rim_c;
-    // As for the 1D case, Vot is not used
-    //double Vot = (length_a * length_b * length_c +
-    //              2.0 * thick_rim_a * length_b * length_c +
-    //              2.0 * length_a * thick_rim_b * length_c +
-    //              2.0 * length_a * length_b * thick_rim_c);
+    const double tA = length_a + 2.0*thick_rim_a;
+    const double tB = length_b + 2.0*thick_rim_b;
+    const double tC = length_c + 2.0*thick_rim_c;
+    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
+    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
+    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
+    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
+    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
+    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
 
-    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
-    // the scaling by B.
-    double ta = length_a + 2.0*thick_rim_a;
-    double tb = length_b + 2.0*thick_rim_b;
-    double tc = length_c + 2.0*thick_rim_c;
-    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    double siA = sas_sinx_x(0.5*length_a*qa);
-    double siB = sas_sinx_x(0.5*length_b*qb);
-    double siC = sas_sinx_x(0.5*length_c*qc);
-    double siAt = sas_sinx_x(0.5*ta*qa);
-    double siBt = sas_sinx_x(0.5*tb*qb);
-    double siCt = sas_sinx_x(0.5*tc*qc);
-
-
-    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
-    // in the 1D code, but should be checked!
-    double f = ( dr0*siA*siB*siC*Vin
-               + drA*(siAt-siA)*siB*siC*V1
-               + drB*siA*(siBt-siB)*siC*V2
-               + drC*siA*siB*(siCt-siC)*V3);
+#if OVERLAPPING
+    const double f = dr0*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siAt*(siBt-siB)*siC
+        + drC*siAt*siBt*(siCt-siC);
+#else
+    const double f = dr0*siA*siB*siC
+        + drA*(siAt-siA)*siB*siC
+        + drB*siA*(siBt-siB)*siC
+        + drC*siA*siB*(siCt-siC);
+#endif
 
     return 1.0e-4 * f * f;

sasmodels/models/core_shell_parallelepiped.py

@@ -5 +5 @@
 Calculates the form factor for a rectangular solid with a core-shell structure.
 The thickness and the scattering length density of the shell or
-"rim" can be different on each (pair) of faces. However at this time the 1D
-calculation does **NOT** actually calculate a c face rim despite the presence
-of the parameter. Some other aspects of the 1D calculation may be wrong.
-
-.. note::
-   This model was originally ported from NIST IGOR macros. However, it is not
-   yet fully understood by the SasView developers and is currently under review.
+"rim" can be different on each (pair) of faces.
 
 The form factor is normalized by the particle volume $V$ such that
@@ -21 +15 @@
 where $\langle \ldots \rangle$ is an average over all possible orientations
 of the rectangular solid.
-
 
 The function calculated is the form factor of the rectangular solid below.
@@ -41 +34 @@
     V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
 
-**meaning that there are "gaps" at the corners of the solid.**  Again note that
-$t_C = 0$ currently.
+**meaning that there are "gaps" at the corners of the solid.**
 
 The intensity calculated follows the :ref:`parallelepiped` model, with the
 core-shell intensity being calculated as the square of the sum of the
-amplitudes of the core and shell, in the same manner as a core-shell model.
-
-.. math::
-
-    F_{a}(Q,\alpha,\beta)=
-    \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)
-               }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta}
-    - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)
-           }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right]
-    \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)
-               }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right]
-    \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)
-               }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right]
-
-.. note::
-
-    Why does t_B not appear in the above equation?
-    For the calculation of the form factor to be valid, the sides of the solid
-    MUST (perhaps not any more?) be chosen such that** $A < B < C$.
-    If this inequality is not satisfied, the model will not report an error,
-    but the calculation will not be correct and thus the result wrong.
+amplitudes of the core and the slabs on the edges.
+
+the scattering amplitude is computed for a particular orientation of the
+core-shell parallelepiped with respect to the scattering vector and then
+averaged over all possible orientations, where $\alpha$ is the angle between
+the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is
+the angle between projection of the particle in the $xy$ detector plane
+and the $y$ axis.
+
+.. math::
+
+    F(Q)
+    &= (\rho_\text{core}-\rho_\text{solvent})
+       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\
+    &+ (\rho_\text{A}-\rho_\text{solvent})
+        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\
+    &+ (\rho_\text{B}-\rho_\text{solvent})
+        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\
+    &+ (\rho_\text{C}-\rho_\text{solvent})
+        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right]
+
+with
+
+.. math::
+
+    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L}
+
+and
+
+.. math::
+
+    Q_A &= \sin\alpha \sin\beta \\
+    Q_B &= \sin\alpha \cos\beta \\
+    Q_C &= \cos\alpha
+
+
+where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$
+are the scattering length of the parallelepiped core, and the rectangular
+slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$
+is the scattering length of the solvent.
 
 FITTING NOTES
+~~~~~~~~~~~~~
+
 If the scale is set equal to the particle volume fraction, $\phi$, the returned
-value is the scattered intensity per unit volume, $I(q) = \phi P(q)$.
-However, **no interparticle interference effects are included in this
-calculation.**
+value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However,
+**no interparticle interference effects are included in this calculation.**
 
 There are many parameters in this model. Hold as many fixed as possible with
 known values, or you will certainly end up at a solution that is unphysical.
 
-Constraints must be applied during fitting to ensure that the inequality
-$A < B < C$ is not violated. The calculation will not report an error,
-but the results will not be correct.
-
 The returned value is in units of |cm^-1|, on absolute scale.
 
 NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated
 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
-and length $(C+2t_C)$ values, after appropriately
-sorting the three dimensions to give an oblate or prolate particle, to give an
-effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
+and length $(C+2t_C)$ values, after appropriately sorting the three dimensions
+to give an oblate or prolate particle, to give an effective radius,
+for $S(Q)$ when $P(Q) * S(Q)$ is applied.
 
 For 2d data the orientation of the particle is required, described using
-angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
-of the calculation and angular dispersions see :ref:`orientation` .
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further
+details of the calculation and angular dispersions see :ref:`orientation`.
 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example,
 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector.
+
+For 2d, constraints must be applied during fitting to ensure that the
+inequality $A < B < C$ is not violated, and hence the correct definition
+of angles is preserved. The calculation will not report an error,
+but the results may be not correct.
 
 .. figure:: img/parallelepiped_angle_definition.png
@@ -114 +127 @@
     Equations (1), (13-14). (in German)
 .. [#] D Singh (2009). *Small angle scattering studies of self assembly in
-   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available
+   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available
    from Proquest <http://search.proquest.com/docview/304915826?accountid
    =26379>`_
@@ -123 +136 @@
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
-* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
-* **Currently Under review by:** Paul Butler
+* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
+* Cross-checked against hollow rectangular prism and rectangular prism for
+  equal thickness overlapping sides, and by Monte Carlo sampling of points
+  within the shape for non-uniform, non-overlapping sides.
 """
 
@@ -175 +190 @@
         Return equivalent radius (ER)
     """
-
-    # surface average radius (rough approximation)
-    surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi)
-
-    height = length_c + 2.0*thick_rim_c
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad))
-    return 0.5 * (ddd) ** (1. / 3.)
+    from .parallelepiped import ER as ER_p
+
+    a = length_a + 2*thick_rim_a
+    b = length_b + 2*thick_rim_b
+    c = length_c + 2*thick_rim_c
+    return ER_p(a, b, c)
 
 # VR defaults to 1.0
@@ -216 +229 @@
             psi_pd=10, psi_pd_n=1)
 
-# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+# rkh 7/4/17 add random unit test for 2d, note make all params different,
+# 2d values not tested against other codes or models
 if 0:  # pak: model rewrite; need to update tests
     qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)

sasmodels/models/cylinder.c

@@ -21 +21 @@
 
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double theta = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N ;i++) {
+        const double theta = GAUSS_Z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         // theta (theta,phi) the projection of the cylinder on the detector plane
         SINCOS(theta , sin_theta, cos_theta);
         const double form = fq(q*sin_theta, q*cos_theta, radius, length);
-        total += Gauss76Wt[i] * form * form * sin_theta;
+        total += GAUSS_W[i] * form * form * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -45 +45 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld,
     double solvent_sld,

sasmodels/models/ellipsoid.c

@@ -22 +22 @@
 
     // translate a point in [-1,1] to a point in [0, 1]
-    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2;
+    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
     const double zm = 0.5;
     const double zb = 0.5;
     double total = 0.0;
-    for (int i=0;i<76;i++) {
-        const double u = Gauss76Z[i]*zm + zb;
+    for (int i=0;i<GAUSS_N;i++) {
+        const double u = GAUSS_Z[i]*zm + zb;
         const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
         const double f = sas_3j1x_x(q*r);
-        total += Gauss76Wt[i] * f * f;
+        total += GAUSS_W[i] * f * f;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -39 +39 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double sld,
     double sld_solvent,

sasmodels/models/elliptical_cylinder.c

@@ -22 +22 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
-        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
         const double sin_val = sqrt(1.0 - cos_val*cos_val);
         //const double arg = radius_minor*sin_val;
         double inner_sum=0;
-        for(int j=0;j<76;j++) {
-            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017
-            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double r = sin_val*sqrt(rA - rB*cos(theta));
             const double be = sas_2J1x_x(q*r);
-            inner_sum += Gauss76Wt[j] * be * be;
+            inner_sum += GAUSS_W[j] * be * be;
         }
         //now calculate the value of the inner integral
@@ -40 +39 @@
         //now calculate outer integral
         const double si = sas_sinx_x(q*0.5*length*cos_val);
-        outer_sum += Gauss76Wt[i] * inner_sum * si * si;
+        outer_sum += GAUSS_W[i] * inner_sum * si * si;
     }
     outer_sum *= 0.5*(vb-va);
@@ -55 +54 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
      double radius_minor, double r_ratio, double length,
      double sld, double solvent_sld)

sasmodels/models/elliptical_cylinder.py

@@ -121 +121 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c",
-          "elliptical_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"]
 
 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0,

sasmodels/models/fcc_paracrystal.c

@@ -53 +53 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -67 +67 @@
             const double qb = qab*sin_phi;
             const double form = fcc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;
@@ -80 +80 @@
 
 
-static double Iqxy(double qa, double qb, double qc,
+static double Iqabc(double qa, double qb, double qc,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld)

sasmodels/models/flexible_cylinder_elliptical.c

@@ -17 +17 @@
     double sum=0.0;
 
-    for(int i=0;i<N_POINTS_76;i++) {
-        const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4;
+    for(int i=0;i<GAUSS_N;i++) {
+        const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4;
         double sn, cn;
         SINCOS(zi, sn, cn);
         const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn);
         const double yyy = sas_2J1x_x(arg);
-        sum += Gauss76Wt[i] * yyy * yyy;
+        sum += GAUSS_W[i] * yyy * yyy;
     }
     sum *= 0.5;

sasmodels/models/hollow_cylinder.c

@@ -38 +38 @@
 
     double summ = 0.0;            //initialize intergral
-    for (int i=0;i<76;i++) {
-        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
+    for (int i=0;i<GAUSS_N;i++) {
+        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper );
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         const double form = _fq(q*sin_theta, q*cos_theta,
                                 radius, thickness, length);
-        summ += Gauss76Wt[i] * form * form;
+        summ += GAUSS_W[i] * form * form;
     }
 
@@ -52 +52 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double radius, double thickness, double length,
     double sld, double solvent_sld)

sasmodels/models/hollow_rectangular_prism.c

@@ -1 +1 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio, double thickness);
 
@@ -37 +37 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
 
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
@@ -49 +49 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
+        for(int j=0; j<GAUSS_N; j++) {
 
-            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -66 +66 @@
             const double AP2 = vol_core * termA2 * termB2 * termC2;
 
-            inner_sum += Gauss76Wt[j] * square(AP1-AP2);
+            inner_sum += GAUSS_W[j] * square(AP1-AP2);
         }
         inner_sum *= 0.5 * (v2b-v2a);
 
-        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta);
+        outer_sum += GAUSS_W[i] * inner_sum * sin(theta);
     }
     outer_sum *= 0.5*(v1b-v1a);
@@ -84 +84 @@
     return 1.0e-4 * delrho * delrho * form;
 }
+
+double Iqabc(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double b2a_ratio,
+    double c2a_ratio,
+    double thickness)
+{
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
+    const double vol_total = length_a * length_b * length_c;
+    const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness);
+
+    // Amplitude AP from eqn. (13)
+
+    const double termA1 = sas_sinx_x(qa * a_half);
+    const double termA2 = sas_sinx_x(qa * (a_half-thickness));
+
+    const double termB1 = sas_sinx_x(qb * b_half);
+    const double termB2 = sas_sinx_x(qb * (b_half-thickness));
+
+    const double termC1 = sas_sinx_x(qc * c_half);
+    const double termC2 = sas_sinx_x(qc * (c_half-thickness));
+
+    const double AP1 = vol_total * termA1 * termB1 * termC1;
+    const double AP2 = vol_core * termA2 * termB2 * termC2;
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    const double delrho = sld - solvent_sld;
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * square(delrho * (AP1-AP2));
+}

sasmodels/models/hollow_rectangular_prism.py

@@ -5 +5 @@
 This model provides the form factor, $P(q)$, for a hollow rectangular
 parallelepiped with a wall of thickness $\Delta$.
-It computes only the 1D scattering, not the 2D.
+
 
 Definition
@@ -66 +66 @@
 (which is unitless).
 
-**The 2D scattering intensity is not computed by this model.**
+For 2d data the orientation of the particle is required, described using
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
+of the calculation and angular dispersions see :ref:`orientation` .
+The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
+$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
+
+For 2d, constraints must be applied during fitting to ensure that the inequality
+$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
+but the results may be not correct.
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented hollow rectangular prism.
+    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
+    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism.
+    The neutron or X-ray beam is along the $z$ axis.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for oriented hollow rectangular prisms against the
+    detector plane.
 
 
@@ -113 +133 @@
               ["thickness", "Ang", 1, [0, inf], "volume",
                "Thickness of parallelepiped"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 

sasmodels/models/hollow_rectangular_prism_thin_walls.c

@@ -1 +1 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio);
 
@@ -29 +29 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
 
         double sin_theta, cos_theta;
@@ -44 +44 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
 
             double sin_phi, cos_phi;
@@ -62 +62 @@
                 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi );
 
-            inner_sum += Gauss76Wt[j] * square(AL+AT);
+            inner_sum += GAUSS_W[j] * square(AL+AT);
         }
 
         inner_sum *= 0.5 * (v2b-v2a);
-        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
 

sasmodels/models/lib/gauss150.c

@@ -1 @@
-/*
- *  GaussWeights.c
- *  SANSAnalysis
- *
- *  Created by Andrew Jackson on 4/23/07.
- *  Copyright 2007 __MyCompanyName__. All rights reserved.
- *
- */
+// Created by Andrew Jackson on 4/23/07
 
-// Gaussians
-constant double Gauss150Z[150]={
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 150
+ #define GAUSS_Z Gauss150Z
+ #define GAUSS_W Gauss150Wt
+
+
+// Note: using array size 152 rather than 150 so that it is a multiple of 4.
+// Some OpenCL devices prefer that vectors start and end on nice boundaries.
+constant double Gauss150Z[152]={
         -0.9998723404457334,
         -0.9993274305065947,
@@ -159 +163 @@
         0.9983473449340834,
         0.9993274305065947,
-        0.9998723404457334
+        0.9998723404457334,
+        0., // zero padding is ignored
+        0.  // zero padding is ignored
 };
 
-constant double Gauss150Wt[150]={
+constant double Gauss150Wt[152]={
         0.0003276086705538,
         0.0007624720924706,
@@ -312 +318 @@
         0.0011976474864367,
         0.0007624720924706,
-        0.0003276086705538
+        0.0003276086705538,
+        0., // zero padding is ignored
+        0.  // zero padding is ignored
 };

sasmodels/models/lib/gauss20.c

@@ -1 @@
-/*
- *  GaussWeights.c
- *  SANSAnalysis
- *
- *  Created by Andrew Jackson on 4/23/07.
- *  Copyright 2007 __MyCompanyName__. All rights reserved.
- *
- */
+// Created by Andrew Jackson on 4/23/07
+
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 20
+ #define GAUSS_Z Gauss20Z
+ #define GAUSS_W Gauss20Wt
 
 // Gaussians

sasmodels/models/lib/gauss76.c

@@ -1 @@
-/*
- *  GaussWeights.c
- *  SANSAnalysis
- *
- *  Created by Andrew Jackson on 4/23/07.
- *  Copyright 2007 __MyCompanyName__. All rights reserved.
- *
- */
-#define N_POINTS_76 76
+// Created by Andrew Jackson on 4/23/07
+
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 76
+ #define GAUSS_Z Gauss76Z
+ #define GAUSS_W Gauss76Wt
 
 // Gaussians
-constant double Gauss76Wt[N_POINTS_76]={
+constant double Gauss76Wt[76]={
         .00126779163408536,             //0
         .00294910295364247,
@@ -89 +90 @@
 };
 
-constant double Gauss76Z[N_POINTS_76]={
+constant double Gauss76Z[76]={
         -.999505948362153,              //0
         -.997397786355355,

sasmodels/models/line.py

@@ -57 +57 @@
 Iq.vectorized = True # Iq accepts an array of q values
 
+
 def Iqxy(qx, qy, *args):
     """
@@ -69 +70 @@
 
 Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values
+
+# uncomment the following to test Iqxy in C models
+#del Iq, Iqxy
+#c_code = """
+#static double Iq(double q, double b, double m) { return m*q+b; }
+#static double Iqxy(double qx, double qy, double b, double m)
+#{ return (m*qx+b)*(m*qy+b); }
+#"""
 
 def random():

sasmodels/models/parallelepiped.c

@@ -23 +23 @@
     double outer_total = 0; //initialize integral
 
-    for( int i=0; i<76; i++) {
-        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
+    for( int i=0; i<GAUSS_N; i++) {
+        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
         const double mu_proj = mu * sqrt(1.0-sigma*sigma);
 
@@ -30 +30 @@
         // corresponding to angles from 0 to pi/2.
         double inner_total = 0.0;
-        for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
             const double si2 = sas_sinx_x(mu_proj * cos_uu);
-            inner_total += Gauss76Wt[j] * square(si1 * si2);
+            inner_total += GAUSS_W[j] * square(si1 * si2);
         }
         inner_total *= 0.5;
 
         const double si = sas_sinx_x(mu * c_scaled * sigma);
-        outer_total += Gauss76Wt[i] * inner_total * si * si;
+        outer_total += GAUSS_W[i] * inner_total * si * si;
     }
     outer_total *= 0.5;
@@ -53 +53 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double sld,
     double solvent_sld,

sasmodels/models/pringle.c

@@ -29 +29 @@
     double sumC = 0.0;          // initialize integral
     double r;
-    for (int i=0; i < 76; i++) {
-        r = Gauss76Z[i]*zm + zb;
+    for (int i=0; i < GAUSS_N; i++) {
+        r = GAUSS_Z[i]*zm + zb;
 
         const double qrs = r*q_sin_psi;
         const double qrrc = r*r*q_cos_psi;
 
-        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
+        double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
         double S, C;
         SINCOS(alpha*qrrc, S, C);
@@ -86 +86 @@
 
     double sum = 0.0;
-    for (int i = 0; i < 76; i++) {
-        double psi = Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++) {
+        double psi = GAUSS_Z[i]*zm + zb;
         double sin_psi, cos_psi;
         SINCOS(psi, sin_psi, cos_psi);
@@ -93 +93 @@
         double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0));
         double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
-        sum += Gauss76Wt[i] * pringle_kernel;
+        sum += GAUSS_W[i] * pringle_kernel;
     }
 

sasmodels/models/rectangular_prism.c

@@ -1 +1 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio);
 
@@ -26 +26 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
@@ -36 +36 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+        for(int j=0; j<GAUSS_N; j++) {
+            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -45 +45 @@
             const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi);
             const double AP = termA * termB * termC;
-            inner_sum += Gauss76Wt[j] * AP * AP;
+            inner_sum += GAUSS_W[j] * AP * AP;
         }
         inner_sum = 0.5 * (v2b-v2a) * inner_sum;
-        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
 
     double answer = 0.5*(v1b-v1a)*outer_sum;
 
-    // Normalize by Pi (Eqn. 16). 
-    // The term (ABC)^2 does not appear because it was introduced before on 
+    // Normalize by Pi (Eqn. 16).
+    // The term (ABC)^2 does not appear because it was introduced before on
     // the definitions of termA, termB, termC.
-    // The factor 2 appears because the theta integral has been defined between 
+    // The factor 2 appears because the theta integral has been defined between
     // 0 and pi/2, instead of 0 to pi.
     answer /= M_PI_2; //Form factor P(q)
@@ -64 +64 @@
     answer *= square((sld-solvent_sld)*volume);
 
-    // Convert from [1e-12 A-1] to [cm-1] 
+    // Convert from [1e-12 A-1] to [cm-1]
     answer *= 1.0e-4;
 
     return answer;
 }
+
+
+double Iqabc(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
+    const double volume = length_a * length_b * length_c;
+
+    // Amplitude AP from eqn. (13)
+
+    const double termA = sas_sinx_x(qa * a_half);
+    const double termB = sas_sinx_x(qb * b_half);
+    const double termC = sas_sinx_x(qc * c_half);
+
+    const double AP = termA * termB * termC;
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    const double delrho = sld - solvent_sld;
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * square(volume * delrho * AP);
+}

sasmodels/models/rectangular_prism.py

@@ -12 +12 @@
 the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity
 to *a* will generate a distribution of cubes of different sizes).
-Note also that, contrary to :ref:`parallelepiped`, it does not compute
-the 2D scattering.
 
 
@@ -26 +24 @@
 that reference), with $\theta$ corresponding to $\alpha$ in that paper,
 and not to the usual convention used for example in the
-:ref:`parallelepiped` model. As the present model does not compute
-the 2D scattering, this has no further consequences.
+:ref:`parallelepiped` model.
 
 In this model the scattering from a massive parallelepiped with an
@@ -65 +62 @@
 units) *scale* represents the volume fraction (which is unitless).
 
-**The 2D scattering intensity is not computed by this model.**
+For 2d data the orientation of the particle is required, described using
+angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details
+of the calculation and angular dispersions see :ref:`orientation` .
+The angle $\Psi$ is the rotational angle around the long *C* axis. For example,
+$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector.
+
+For 2d, constraints must be applied during fitting to ensure that the inequality
+$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error,
+but the results may be not correct.
+
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented core-shell parallelepipeds.
+    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then
+    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
+    The neutron or X-ray beam is along the $z$ axis.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for oriented rectangular prisms against the
+    detector plane.
+
 
 
@@ -108 +126 @@
               ["c2a_ratio", "", 1, [0, inf], "volume",
                "Ratio sides c/a"],
+              ["theta", "degrees", 0, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 0, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 

sasmodels/models/sc_paracrystal.c

@@ -54 +54 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -68 +68 @@
             const double qb = qab*sin_phi;
             const double form = sc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;
@@ -82 +82 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double dnn, double d_factor, double radius,
     double sld, double solvent_sld)

sasmodels/models/stacked_disks.c

@@ -81 +81 @@
     double halfheight = 0.5*thick_core;
 
-    for(int i=0; i<N_POINTS_76; i++) {
-        double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
+    for(int i=0; i<GAUSS_N; i++) {
+        double zi = (GAUSS_Z[i] + 1.0)*M_PI_4;
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(zi, sin_alpha, cos_alpha);
@@ -95 +95 @@
                            solvent_sld,
                            d);
-        summ += Gauss76Wt[i] * yyy * sin_alpha;
+        summ += GAUSS_W[i] * yyy * sin_alpha;
     }
 
@@ -142 +142 @@
 
 static double
-Iqxy(double qab, double qc,
+Iqac(double qab, double qc,
     double thick_core,
     double thick_layer,

sasmodels/models/triaxial_ellipsoid.c

@@ -21 +21 @@
     const double zb = M_PI_4;
     double outer = 0.0;
-    for (int i=0;i<76;i++) {
-        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2;
-        const double phi = Gauss76Z[i]*zm + zb;
+    for (int i=0;i<GAUSS_N;i++) {
+        //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2;
+        const double phi = GAUSS_Z[i]*zm + zb;
         const double pa_sinsq_phi = pa*square(sin(phi));
 
@@ -29 +29 @@
         const double um = 0.5;
         const double ub = 0.5;
-        for (int j=0;j<76;j++) {
+        for (int j=0;j<GAUSS_N;j++) {
             // translate a point in [-1,1] to a point in [0, 1]
-            const double usq = square(Gauss76Z[j]*um + ub);
+            const double usq = square(GAUSS_Z[j]*um + ub);
             const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq);
             const double fq = sas_3j1x_x(q*r);
-            inner += Gauss76Wt[j] * fq * fq;
+            inner += GAUSS_W[j] * fq * fq;
         }
-        outer += Gauss76Wt[i] * inner;  // correcting for dx later
+        outer += GAUSS_W[i] * inner;  // correcting for dx later
     }
     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
@@ -46 +46 @@
 
 static double
-Iqxy(double qa, double qb, double qc,
+Iqabc(double qa, double qb, double qc,
     double sld,
     double sld_solvent,

sasmodels/special.py

@@ -3 +3 @@
 .................
 
-The C code follows the C99 standard, with the usual math functions,
-as defined in
-`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_.
-This includes the following:
+This following standard C99 math functions are available:
 
     M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:
         $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$
 
-    exp, log, pow(x,y), expm1, sqrt:
-        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$.
-        The function expm1(x) is accurate across all $x$, including $x$
-        very close to zero.
+    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt:
+        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$,
+        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x)
+        are accurate across all $x$, including $x$ very close to zero.
 
     sin, cos, tan, asin, acos, atan:
@@ -29 +26 @@
         quadrants II and IV when $x$ and $y$ have opposite sign.
 
-    fmin(x,y), fmax(x,y), trunc, rint:
+    fabs(x), fmin(x,y), fmax(x,y), trunc, rint:
         Floating point functions.  rint(x) returns the nearest integer.
 
@@ -35 +32 @@
         NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that
         you cannot use :code:`x == NAN` to test for NaN values since that
-        will always return false.  NAN does not equal NAN!
+        will always return false.  NAN does not equal NAN!  The alternative,
+        :code:`x != x` may fail if the compiler optimizes the test away.
 
     INFINITY:
@@ -89 +87 @@
         for forcing a constant to stay double precision.
 
-The following special functions and scattering calculations are defined in
-`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_.
+The following special functions and scattering calculations are defined.
 These functions have been tuned to be fast and numerically stable down
 to $q=0$ even in single precision.  In some cases they work around bugs
@@ -184 +181 @@
 
 
-    Gauss76Z[i], Gauss76Wt[i]:
+    gauss76.n, gauss76.z[i], gauss76.w[i]:
         Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively,
         computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$.
-
-        Similar arrays are available in :code:`gauss20.c` for 20-point
-        quadrature and in :code:`gauss150.c` for 150-point quadrature.
-
+        When translating the model to C, include 'lib/gauss76.c' in the source
+        and use :code:`GAUSS_N`, :code:`GAUSS_Z`, and :code:`GAUSS_W`.
+
+        Similar arrays are available in :code:`gauss20` for 20-point quadrature
+        and :code:`gauss150.c` for 150-point quadrature. By using
+        :code:`import gauss76 as gauss` it is easy to change the number of
+        points in the integration.
 """
 # pylint: disable=unused-import
@@ -200 +200 @@
 
 # C99 standard math library functions
-from numpy import exp, log, power as pow, expm1, sqrt
+from numpy import exp, log, power as pow, expm1, log1p, sqrt, cbrt
 from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan
 from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh
 from numpy import arctan2 as atan2
-from numpy import fmin, fmax, trunc, rint
-from numpy import NAN, inf as INFINITY
-
+from numpy import fabs, fmin, fmax, trunc, rint
+from numpy import pi, nan, inf
 from scipy.special import gamma as sas_gamma
 from scipy.special import erf as sas_erf
@@ -218 +217 @@
 # C99 standard math constants
 M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E = np.pi, np.pi/2, np.pi/4, np.sqrt(0.5), np.e
+NAN = nan
+INFINITY = inf
 
 # non-standard constants
@@ -226 +227 @@
     """return sin(x), cos(x)"""
     return sin(x), cos(x)
+sincos = SINCOS
 
 def square(x):
@@ -294 +296 @@
 
 # Gaussians
-
-Gauss20Wt = np.array([
-    .0176140071391521,
-    .0406014298003869,
-    .0626720483341091,
-    .0832767415767047,
-    .10193011981724,
-    .118194531961518,
-    .131688638449177,
-    .142096109318382,
-    .149172986472604,
-    .152753387130726,
-    .152753387130726,
-    .149172986472604,
-    .142096109318382,
-    .131688638449177,
-    .118194531961518,
-    .10193011981724,
-    .0832767415767047,
-    .0626720483341091,
-    .0406014298003869,
-    .0176140071391521
-])
-
-Gauss20Z = np.array([
-    -.993128599185095,
-    -.963971927277914,
-    -.912234428251326,
-    -.839116971822219,
-    -.746331906460151,
-    -.636053680726515,
-    -.510867001950827,
-    -.37370608871542,
-    -.227785851141645,
-    -.076526521133497,
-    .0765265211334973,
-    .227785851141645,
-    .37370608871542,
-    .510867001950827,
-    .636053680726515,
-    .746331906460151,
-    .839116971822219,
-    .912234428251326,
-    .963971927277914,
-    .993128599185095
-])
-
-Gauss76Wt = np.array([
-    .00126779163408536,         #0
-    .00294910295364247,
-    .00462793522803742,
-    .00629918049732845,
-    .00795984747723973,
-    .00960710541471375,
-    .0112381685696677,
-    .0128502838475101,
-    .0144407317482767,
-    .0160068299122486,
-    .0175459372914742,          #10
-    .0190554584671906,
-    .020532847967908,
-    .0219756145344162,
-    .0233813253070112,
-    .0247476099206597,
-    .026072164497986,
-    .0273527555318275,
-    .028587223650054,
-    .029773487255905,
-    .0309095460374916,          #20
-    .0319934843404216,
-    .0330234743977917,
-    .0339977794120564,
-    .0349147564835508,
-    .0357728593807139,
-    .0365706411473296,
-    .0373067565423816,
-    .0379799643084053,
-    .0385891292645067,
-    .0391332242205184,          #30
-    .0396113317090621,
-    .0400226455325968,
-    .040366472122844,
-    .0406422317102947,
-    .0408494593018285,
-    .040987805464794,
-    .0410570369162294,
-    .0410570369162294,
-    .040987805464794,
-    .0408494593018285,          #40
-    .0406422317102947,
-    .040366472122844,
-    .0400226455325968,
-    .0396113317090621,
-    .0391332242205184,
-    .0385891292645067,
-    .0379799643084053,
-    .0373067565423816,
-    .0365706411473296,
-    .0357728593807139,          #50
-    .0349147564835508,
-    .0339977794120564,
-    .0330234743977917,
-    .0319934843404216,
-    .0309095460374916,
-    .029773487255905,
-    .028587223650054,
-    .0273527555318275,
-    .026072164497986,
-    .0247476099206597,          #60
-    .0233813253070112,
-    .0219756145344162,
-    .020532847967908,
-    .0190554584671906,
-    .0175459372914742,
-    .0160068299122486,
-    .0144407317482767,
-    .0128502838475101,
-    .0112381685696677,
-    .00960710541471375,         #70
-    .00795984747723973,
-    .00629918049732845,
-    .00462793522803742,
-    .00294910295364247,
-    .00126779163408536          #75 (indexed from 0)
-])
-
-Gauss76Z = np.array([
-    -.999505948362153,          #0
-    -.997397786355355,
-    -.993608772723527,
-    -.988144453359837,
-    -.981013938975656,
-    -.972229228520377,
-    -.961805126758768,
-    -.949759207710896,
-    -.936111781934811,
-    -.92088586125215,
-    -.904107119545567,          #10
-    -.885803849292083,
-    -.866006913771982,
-    -.844749694983342,
-    -.822068037328975,
-    -.7980001871612,
-    -.77258672828181,
-    -.74587051350361,
-    -.717896592387704,
-    -.688712135277641,
-    -.658366353758143,          #20
-    -.626910417672267,
-    -.594397368836793,
-    -.560882031601237,
-    -.526420920401243,
-    -.491072144462194,
-    -.454895309813726,
-    -.417951418780327,
-    -.380302767117504,
-    -.342012838966962,
-    -.303146199807908,          #30
-    -.263768387584994,
-    -.223945802196474,
-    -.183745593528914,
-    -.143235548227268,
-    -.102483975391227,
-    -.0615595913906112,
-    -.0205314039939986,
-    .0205314039939986,
-    .0615595913906112,
-    .102483975391227,                   #40
-    .143235548227268,
-    .183745593528914,
-    .223945802196474,
-    .263768387584994,
-    .303146199807908,
-    .342012838966962,
-    .380302767117504,
-    .417951418780327,
-    .454895309813726,
-    .491072144462194,           #50
-    .526420920401243,
-    .560882031601237,
-    .594397368836793,
-    .626910417672267,
-    .658366353758143,
-    .688712135277641,
-    .717896592387704,
-    .74587051350361,
-    .77258672828181,
-    .7980001871612,     #60
-    .822068037328975,
-    .844749694983342,
-    .866006913771982,
-    .885803849292083,
-    .904107119545567,
-    .92088586125215,
-    .936111781934811,
-    .949759207710896,
-    .961805126758768,
-    .972229228520377,           #70
-    .981013938975656,
-    .988144453359837,
-    .993608772723527,
-    .997397786355355,
-    .999505948362153            #75
-])
-
-Gauss150Z = np.array([
-    -0.9998723404457334,
-    -0.9993274305065947,
-    -0.9983473449340834,
-    -0.9969322929775997,
-    -0.9950828645255290,
-    -0.9927998590434373,
-    -0.9900842691660192,
-    -0.9869372772712794,
-    -0.9833602541697529,
-    -0.9793547582425894,
-    -0.9749225346595943,
-    -0.9700655145738374,
-    -0.9647858142586956,
-    -0.9590857341746905,
-    -0.9529677579610971,
-    -0.9464345513503147,
-    -0.9394889610042837,
-    -0.9321340132728527,
-    -0.9243729128743136,
-    -0.9162090414984952,
-    -0.9076459563329236,
-    -0.8986873885126239,
-    -0.8893372414942055,
-    -0.8795995893549102,
-    -0.8694786750173527,
-    -0.8589789084007133,
-    -0.8481048644991847,
-    -0.8368612813885015,
-    -0.8252530581614230,
-    -0.8132852527930605,
-    -0.8009630799369827,
-    -0.7882919086530552,
-    -0.7752772600680049,
-    -0.7619248049697269,
-    -0.7482403613363824,
-    -0.7342298918013638,
-    -0.7198995010552305,
-    -0.7052554331857488,
-    -0.6903040689571928,
-    -0.6750519230300931,
-    -0.6595056411226444,
-    -0.6436719971150083,
-    -0.6275578900977726,
-    -0.6111703413658551,
-    -0.5945164913591590,
-    -0.5776035965513142,
-    -0.5604390262878617,
-    -0.5430302595752546,
-    -0.5253848818220803,
-    -0.5075105815339176,
-    -0.4894151469632753,
-    -0.4711064627160663,
-    -0.4525925063160997,
-    -0.4338813447290861,
-    -0.4149811308476706,
-    -0.3959000999390257,
-    -0.3766465660565522,
-    -0.3572289184172501,
-    -0.3376556177463400,
-    -0.3179351925907259,
-    -0.2980762356029071,
-    -0.2780873997969574,
-    -0.2579773947782034,
-    -0.2377549829482451,
-    -0.2174289756869712,
-    -0.1970082295132342,
-    -0.1765016422258567,
-    -0.1559181490266516,
-    -0.1352667186271445,
-    -0.1145563493406956,
-    -0.0937960651617229,
-    -0.0729949118337358,
-    -0.0521619529078925,
-    -0.0313062657937972,
-    -0.0104369378042598,
-    0.0104369378042598,
-    0.0313062657937972,
-    0.0521619529078925,
-    0.0729949118337358,
-    0.0937960651617229,
-    0.1145563493406956,
-    0.1352667186271445,
-    0.1559181490266516,
-    0.1765016422258567,
-    0.1970082295132342,
-    0.2174289756869712,
-    0.2377549829482451,
-    0.2579773947782034,
-    0.2780873997969574,
-    0.2980762356029071,
-    0.3179351925907259,
-    0.3376556177463400,
-    0.3572289184172501,
-    0.3766465660565522,
-    0.3959000999390257,
-    0.4149811308476706,
-    0.4338813447290861,
-    0.4525925063160997,
-    0.4711064627160663,
-    0.4894151469632753,
-    0.5075105815339176,
-    0.5253848818220803,
-    0.5430302595752546,
-    0.5604390262878617,
-    0.5776035965513142,
-    0.5945164913591590,
-    0.6111703413658551,
-    0.6275578900977726,
-    0.6436719971150083,
-    0.6595056411226444,
-    0.6750519230300931,
-    0.6903040689571928,
-    0.7052554331857488,
-    0.7198995010552305,
-    0.7342298918013638,
-    0.7482403613363824,
-    0.7619248049697269,
-    0.7752772600680049,
-    0.7882919086530552,
-    0.8009630799369827,
-    0.8132852527930605,
-    0.8252530581614230,
-    0.8368612813885015,
-    0.8481048644991847,
-    0.8589789084007133,
-    0.8694786750173527,
-    0.8795995893549102,
-    0.8893372414942055,
-    0.8986873885126239,
-    0.9076459563329236,
-    0.9162090414984952,
-    0.9243729128743136,
-    0.9321340132728527,
-    0.9394889610042837,
-    0.9464345513503147,
-    0.9529677579610971,
-    0.9590857341746905,
-    0.9647858142586956,
-    0.9700655145738374,
-    0.9749225346595943,
-    0.9793547582425894,
-    0.9833602541697529,
-    0.9869372772712794,
-    0.9900842691660192,
-    0.9927998590434373,
-    0.9950828645255290,
-    0.9969322929775997,
-    0.9983473449340834,
-    0.9993274305065947,
-    0.9998723404457334
-])
-
-Gauss150Wt = np.array([
-    0.0003276086705538,
-    0.0007624720924706,
-    0.0011976474864367,
-    0.0016323569986067,
-    0.0020663664924131,
-    0.0024994789888943,
-    0.0029315036836558,
-    0.0033622516236779,
-    0.0037915348363451,
-    0.0042191661429919,
-    0.0046449591497966,
-    0.0050687282939456,
-    0.0054902889094487,
-    0.0059094573005900,
-    0.0063260508184704,
-    0.0067398879387430,
-    0.0071507883396855,
-    0.0075585729801782,
-    0.0079630641773633,
-    0.0083640856838475,
-    0.0087614627643580,
-    0.0091550222717888,
-    0.0095445927225849,
-    0.0099300043714212,
-    0.0103110892851360,
-    0.0106876814158841,
-    0.0110596166734735,
-    0.0114267329968529,
-    0.0117888704247183,
-    0.0121458711652067,
-    0.0124975796646449,
-    0.0128438426753249,
-    0.0131845093222756,
-    0.0135194311690004,
-    0.0138484622795371,
-    0.0141714592928592,
-    0.0144882814685445,
-    0.0147987907597169,
-    0.0151028518701744,
-    0.0154003323133401,
-    0.0156911024699895,
-    0.0159750356447283,
-    0.0162520081211971,
-    0.0165218992159766,
-    0.0167845913311726,
-    0.0170399700056559,
-    0.0172879239649355,
-    0.0175283451696437,
-    0.0177611288626114,
-    0.0179861736145128,
-    0.0182033813680609,
-    0.0184126574807331,
-    0.0186139107660094,
-    0.0188070535331042,
-    0.0189920016251754,
-    0.0191686744559934,
-    0.0193369950450545,
-    0.0194968900511231,
-    0.0196482898041878,
-    0.0197911283358190,
-    0.0199253434079123,
-    0.0200508765398072,
-    0.0201676730337687,
-    0.0202756819988200,
-    0.0203748563729175,
-    0.0204651529434560,
-    0.0205465323660984,
-    0.0206189591819181,
-    0.0206824018328499,
-    0.0207368326754401,
-    0.0207822279928917,
-    0.0208185680053983,
-    0.0208458368787627,
-    0.0208640227312962,
-    0.0208731176389954,
-    0.0208731176389954,
-    0.0208640227312962,
-    0.0208458368787627,
-    0.0208185680053983,
-    0.0207822279928917,
-    0.0207368326754401,
-    0.0206824018328499,
-    0.0206189591819181,
-    0.0205465323660984,
-    0.0204651529434560,
-    0.0203748563729175,
-    0.0202756819988200,
-    0.0201676730337687,
-    0.0200508765398072,
-    0.0199253434079123,
-    0.0197911283358190,
-    0.0196482898041878,
-    0.0194968900511231,
-    0.0193369950450545,
-    0.0191686744559934,
-    0.0189920016251754,
-    0.0188070535331042,
-    0.0186139107660094,
-    0.0184126574807331,
-    0.0182033813680609,
-    0.0179861736145128,
-    0.0177611288626114,
-    0.0175283451696437,
-    0.0172879239649355,
-    0.0170399700056559,
-    0.0167845913311726,
-    0.0165218992159766,
-    0.0162520081211971,
-    0.0159750356447283,
-    0.0156911024699895,
-    0.0154003323133401,
-    0.0151028518701744,
-    0.0147987907597169,
-    0.0144882814685445,
-    0.0141714592928592,
-    0.0138484622795371,
-    0.0135194311690004,
-    0.0131845093222756,
-    0.0128438426753249,
-    0.0124975796646449,
-    0.0121458711652067,
-    0.0117888704247183,
-    0.0114267329968529,
-    0.0110596166734735,
-    0.0106876814158841,
-    0.0103110892851360,
-    0.0099300043714212,
-    0.0095445927225849,
-    0.0091550222717888,
-    0.0087614627643580,
-    0.0083640856838475,
-    0.0079630641773633,
-    0.0075585729801782,
-    0.0071507883396855,
-    0.0067398879387430,
-    0.0063260508184704,
-    0.0059094573005900,
-    0.0054902889094487,
-    0.0050687282939456,
-    0.0046449591497966,
-    0.0042191661429919,
-    0.0037915348363451,
-    0.0033622516236779,
-    0.0029315036836558,
-    0.0024994789888943,
-    0.0020663664924131,
-    0.0016323569986067,
-    0.0011976474864367,
-    0.0007624720924706,
-    0.0003276086705538
-])
+class Gauss:
+    def __init__(self, w, z):
+        self.n = len(w)
+        self.w = w
+        self.z = z
+
+gauss20 = Gauss(
+    w=np.array([
+        .0176140071391521,
+        .0406014298003869,
+        .0626720483341091,
+        .0832767415767047,
+        .10193011981724,
+        .118194531961518,
+        .131688638449177,
+        .142096109318382,
+        .149172986472604,
+        .152753387130726,
+        .152753387130726,
+        .149172986472604,
+        .142096109318382,
+        .131688638449177,
+        .118194531961518,
+        .10193011981724,
+        .0832767415767047,
+        .0626720483341091,
+        .0406014298003869,
+        .0176140071391521
+    ]),
+    z=np.array([
+        -.993128599185095,
+        -.963971927277914,
+        -.912234428251326,
+        -.839116971822219,
+        -.746331906460151,
+        -.636053680726515,
+        -.510867001950827,
+        -.37370608871542,
+        -.227785851141645,
+        -.076526521133497,
+        .0765265211334973,
+        .227785851141645,
+        .37370608871542,
+        .510867001950827,
+        .636053680726515,
+        .746331906460151,
+        .839116971822219,
+        .912234428251326,
+        .963971927277914,
+        .993128599185095
+    ])
+)
+
+gauss76 = Gauss(
+    w=np.array([
+        .00126779163408536,             #0
+        .00294910295364247,
+        .00462793522803742,
+        .00629918049732845,
+        .00795984747723973,
+        .00960710541471375,
+        .0112381685696677,
+        .0128502838475101,
+        .0144407317482767,
+        .0160068299122486,
+        .0175459372914742,              #10
+        .0190554584671906,
+        .020532847967908,
+        .0219756145344162,
+        .0233813253070112,
+        .0247476099206597,
+        .026072164497986,
+        .0273527555318275,
+        .028587223650054,
+        .029773487255905,
+        .0309095460374916,              #20
+        .0319934843404216,
+        .0330234743977917,
+        .0339977794120564,
+        .0349147564835508,
+        .0357728593807139,
+        .0365706411473296,
+        .0373067565423816,
+        .0379799643084053,
+        .0385891292645067,
+        .0391332242205184,              #30
+        .0396113317090621,
+        .0400226455325968,
+        .040366472122844,
+        .0406422317102947,
+        .0408494593018285,
+        .040987805464794,
+        .0410570369162294,
+        .0410570369162294,
+        .040987805464794,
+        .0408494593018285,              #40
+        .0406422317102947,
+        .040366472122844,
+        .0400226455325968,
+        .0396113317090621,
+        .0391332242205184,
+        .0385891292645067,
+        .0379799643084053,
+        .0373067565423816,
+        .0365706411473296,
+        .0357728593807139,              #50
+        .0349147564835508,
+        .0339977794120564,
+        .0330234743977917,
+        .0319934843404216,
+        .0309095460374916,
+        .029773487255905,
+        .028587223650054,
+        .0273527555318275,
+        .026072164497986,
+        .0247476099206597,              #60
+        .0233813253070112,
+        .0219756145344162,
+        .020532847967908,
+        .0190554584671906,
+        .0175459372914742,
+        .0160068299122486,
+        .0144407317482767,
+        .0128502838475101,
+        .0112381685696677,
+        .00960710541471375,             #70
+        .00795984747723973,
+        .00629918049732845,
+        .00462793522803742,
+        .00294910295364247,
+        .00126779163408536              #75 (indexed from 0)
+    ]),
+    z=np.array([
+        -.999505948362153,              #0
+        -.997397786355355,
+        -.993608772723527,
+        -.988144453359837,
+        -.981013938975656,
+        -.972229228520377,
+        -.961805126758768,
+        -.949759207710896,
+        -.936111781934811,
+        -.92088586125215,
+        -.904107119545567,              #10
+        -.885803849292083,
+        -.866006913771982,
+        -.844749694983342,
+        -.822068037328975,
+        -.7980001871612,
+        -.77258672828181,
+        -.74587051350361,
+        -.717896592387704,
+        -.688712135277641,
+        -.658366353758143,              #20
+        -.626910417672267,
+        -.594397368836793,
+        -.560882031601237,
+        -.526420920401243,
+        -.491072144462194,
+        -.454895309813726,
+        -.417951418780327,
+        -.380302767117504,
+        -.342012838966962,
+        -.303146199807908,              #30
+        -.263768387584994,
+        -.223945802196474,
+        -.183745593528914,
+        -.143235548227268,
+        -.102483975391227,
+        -.0615595913906112,
+        -.0205314039939986,
+        .0205314039939986,
+        .0615595913906112,
+        .102483975391227,                       #40
+        .143235548227268,
+        .183745593528914,
+        .223945802196474,
+        .263768387584994,
+        .303146199807908,
+        .342012838966962,
+        .380302767117504,
+        .417951418780327,
+        .454895309813726,
+        .491072144462194,               #50
+        .526420920401243,
+        .560882031601237,
+        .594397368836793,
+        .626910417672267,
+        .658366353758143,
+        .688712135277641,
+        .717896592387704,
+        .74587051350361,
+        .77258672828181,
+        .7980001871612, #60
+        .822068037328975,
+        .844749694983342,
+        .866006913771982,
+        .885803849292083,
+        .904107119545567,
+        .92088586125215,
+        .936111781934811,
+        .949759207710896,
+        .961805126758768,
+        .972229228520377,               #70
+        .981013938975656,
+        .988144453359837,
+        .993608772723527,
+        .997397786355355,
+        .999505948362153                #75
+    ])
+)
+
+gauss150 = Gauss(
+    z=np.array([
+        -0.9998723404457334,
+        -0.9993274305065947,
+        -0.9983473449340834,
+        -0.9969322929775997,
+        -0.9950828645255290,
+        -0.9927998590434373,
+        -0.9900842691660192,
+        -0.9869372772712794,
+        -0.9833602541697529,
+        -0.9793547582425894,
+        -0.9749225346595943,
+        -0.9700655145738374,
+        -0.9647858142586956,
+        -0.9590857341746905,
+        -0.9529677579610971,
+        -0.9464345513503147,
+        -0.9394889610042837,
+        -0.9321340132728527,
+        -0.9243729128743136,
+        -0.9162090414984952,
+        -0.9076459563329236,
+        -0.8986873885126239,
+        -0.8893372414942055,
+        -0.8795995893549102,
+        -0.8694786750173527,
+        -0.8589789084007133,
+        -0.8481048644991847,
+        -0.8368612813885015,
+        -0.8252530581614230,
+        -0.8132852527930605,
+        -0.8009630799369827,
+        -0.7882919086530552,
+        -0.7752772600680049,
+        -0.7619248049697269,
+        -0.7482403613363824,
+        -0.7342298918013638,
+        -0.7198995010552305,
+        -0.7052554331857488,
+        -0.6903040689571928,
+        -0.6750519230300931,
+        -0.6595056411226444,
+        -0.6436719971150083,
+        -0.6275578900977726,
+        -0.6111703413658551,
+        -0.5945164913591590,
+        -0.5776035965513142,
+        -0.5604390262878617,
+        -0.5430302595752546,
+        -0.5253848818220803,
+        -0.5075105815339176,
+        -0.4894151469632753,
+        -0.4711064627160663,
+        -0.4525925063160997,
+        -0.4338813447290861,
+        -0.4149811308476706,
+        -0.3959000999390257,
+        -0.3766465660565522,
+        -0.3572289184172501,
+        -0.3376556177463400,
+        -0.3179351925907259,
+        -0.2980762356029071,
+        -0.2780873997969574,
+        -0.2579773947782034,
+        -0.2377549829482451,
+        -0.2174289756869712,
+        -0.1970082295132342,
+        -0.1765016422258567,
+        -0.1559181490266516,
+        -0.1352667186271445,
+        -0.1145563493406956,
+        -0.0937960651617229,
+        -0.0729949118337358,
+        -0.0521619529078925,
+        -0.0313062657937972,
+        -0.0104369378042598,
+        0.0104369378042598,
+        0.0313062657937972,
+        0.0521619529078925,
+        0.0729949118337358,
+        0.0937960651617229,
+        0.1145563493406956,
+        0.1352667186271445,
+        0.1559181490266516,
+        0.1765016422258567,
+        0.1970082295132342,
+        0.2174289756869712,
+        0.2377549829482451,
+        0.2579773947782034,
+        0.2780873997969574,
+        0.2980762356029071,
+        0.3179351925907259,
+        0.3376556177463400,
+        0.3572289184172501,
+        0.3766465660565522,
+        0.3959000999390257,
+        0.4149811308476706,
+        0.4338813447290861,
+        0.4525925063160997,
+        0.4711064627160663,
+        0.4894151469632753,
+        0.5075105815339176,
+        0.5253848818220803,
+        0.5430302595752546,
+        0.5604390262878617,
+        0.5776035965513142,
+        0.5945164913591590,
+        0.6111703413658551,
+        0.6275578900977726,
+        0.6436719971150083,
+        0.6595056411226444,
+        0.6750519230300931,
+        0.6903040689571928,
+        0.7052554331857488,
+        0.7198995010552305,
+        0.7342298918013638,
+        0.7482403613363824,
+        0.7619248049697269,
+        0.7752772600680049,
+        0.7882919086530552,
+        0.8009630799369827,
+        0.8132852527930605,
+        0.8252530581614230,
+        0.8368612813885015,
+        0.8481048644991847,
+        0.8589789084007133,
+        0.8694786750173527,
+        0.8795995893549102,
+        0.8893372414942055,
+        0.8986873885126239,
+        0.9076459563329236,
+        0.9162090414984952,
+        0.9243729128743136,
+        0.9321340132728527,
+        0.9394889610042837,
+        0.9464345513503147,
+        0.9529677579610971,
+        0.9590857341746905,
+        0.9647858142586956,
+        0.9700655145738374,
+        0.9749225346595943,
+        0.9793547582425894,
+        0.9833602541697529,
+        0.9869372772712794,
+        0.9900842691660192,
+        0.9927998590434373,
+        0.9950828645255290,
+        0.9969322929775997,
+        0.9983473449340834,
+        0.9993274305065947,
+        0.9998723404457334
+    ]),
+    w=np.array([
+        0.0003276086705538,
+        0.0007624720924706,
+        0.0011976474864367,
+        0.0016323569986067,
+        0.0020663664924131,
+        0.0024994789888943,
+        0.0029315036836558,
+        0.0033622516236779,
+        0.0037915348363451,
+        0.0042191661429919,
+        0.0046449591497966,
+        0.0050687282939456,
+        0.0054902889094487,
+        0.0059094573005900,
+        0.0063260508184704,
+        0.0067398879387430,
+        0.0071507883396855,
+        0.0075585729801782,
+        0.0079630641773633,
+        0.0083640856838475,
+        0.0087614627643580,
+        0.0091550222717888,
+        0.0095445927225849,
+        0.0099300043714212,
+        0.0103110892851360,
+        0.0106876814158841,
+        0.0110596166734735,
+        0.0114267329968529,
+        0.0117888704247183,
+        0.0121458711652067,
+        0.0124975796646449,
+        0.0128438426753249,
+        0.0131845093222756,
+        0.0135194311690004,
+        0.0138484622795371,
+        0.0141714592928592,
+        0.0144882814685445,
+        0.0147987907597169,
+        0.0151028518701744,
+        0.0154003323133401,
+        0.0156911024699895,
+        0.0159750356447283,
+        0.0162520081211971,
+        0.0165218992159766,
+        0.0167845913311726,
+        0.0170399700056559,
+        0.0172879239649355,
+        0.0175283451696437,
+        0.0177611288626114,
+        0.0179861736145128,
+        0.0182033813680609,
+        0.0184126574807331,
+        0.0186139107660094,
+        0.0188070535331042,
+        0.0189920016251754,
+        0.0191686744559934,
+        0.0193369950450545,
+        0.0194968900511231,
+        0.0196482898041878,
+        0.0197911283358190,
+        0.0199253434079123,
+        0.0200508765398072,
+        0.0201676730337687,
+        0.0202756819988200,
+        0.0203748563729175,
+        0.0204651529434560,
+        0.0205465323660984,
+        0.0206189591819181,
+        0.0206824018328499,
+        0.0207368326754401,
+        0.0207822279928917,
+        0.0208185680053983,
+        0.0208458368787627,
+        0.0208640227312962,
+        0.0208731176389954,
+        0.0208731176389954,
+        0.0208640227312962,
+        0.0208458368787627,
+        0.0208185680053983,
+        0.0207822279928917,
+        0.0207368326754401,
+        0.0206824018328499,
+        0.0206189591819181,
+        0.0205465323660984,
+        0.0204651529434560,
+        0.0203748563729175,
+        0.0202756819988200,
+        0.0201676730337687,
+        0.0200508765398072,
+        0.0199253434079123,
+        0.0197911283358190,
+        0.0196482898041878,
+        0.0194968900511231,
+        0.0193369950450545,
+        0.0191686744559934,
+        0.0189920016251754,
+        0.0188070535331042,
+        0.0186139107660094,
+        0.0184126574807331,
+        0.0182033813680609,
+        0.0179861736145128,
+        0.0177611288626114,
+        0.0175283451696437,
+        0.0172879239649355,
+        0.0170399700056559,
+        0.0167845913311726,
+        0.0165218992159766,
+        0.0162520081211971,
+        0.0159750356447283,
+        0.0156911024699895,
+        0.0154003323133401,
+        0.0151028518701744,
+        0.0147987907597169,
+        0.0144882814685445,
+        0.0141714592928592,
+        0.0138484622795371,
+        0.0135194311690004,
+        0.0131845093222756,
+        0.0128438426753249,
+        0.0124975796646449,
+        0.0121458711652067,
+        0.0117888704247183,
+        0.0114267329968529,
+        0.0110596166734735,
+        0.0106876814158841,
+        0.0103110892851360,
+        0.0099300043714212,
+        0.0095445927225849,
+        0.0091550222717888,
+        0.0087614627643580,
+        0.0083640856838475,
+        0.0079630641773633,
+        0.0075585729801782,
+        0.0071507883396855,
+        0.0067398879387430,
+        0.0063260508184704,
+        0.0059094573005900,
+        0.0054902889094487,
+        0.0050687282939456,
+        0.0046449591497966,
+        0.0042191661429919,
+        0.0037915348363451,
+        0.0033622516236779,
+        0.0029315036836558,
+        0.0024994789888943,
+        0.0020663664924131,
+        0.0016323569986067,
+        0.0011976474864367,
+        0.0007624720924706,
+        0.0003276086705538
+    ])
+)