# Source code for fluids.particle_size_distribution

    """Chemical Engineering Design Library (ChEDL). Utilities for process modeling.
    Copyright (C) 2016, 2017, 2018, 2019, 2020 Caleb Bell <Caleb.Andrew.Bell@gmail.com>

    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:

    The above copyright notice and this permission notice shall be included in all
    copies or substantial portions of the Software.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
    SOFTWARE.

    This module contains particle distribution characterization, fitting,
    interpolating, and manipulation functions. It may be used with discrete
    particle size distributions, or with statistical ones with parameters
    specified.

    For reporting bugs, adding feature requests, or submitting pull requests,
    please use the `GitHub issue tracker <https://github.com/CalebBell/fluids/>`_
    or contact the author at Caleb.Andrew.Bell@gmail.com.

    .. contents:: :local:

    Particle Size Distribution Base Class
    -------------------------------------
    .. autoclass:: ParticleSizeDistributionContinuous
        :members:

    Discrete Particle Size Distributions
    ------------------------------------
    .. autoclass:: ParticleSizeDistribution
        :members:
    .. autoclass:: PSDInterpolated
        :members:

    Statistical Particle Size Distributions
    ---------------------------------------
    .. autoclass:: PSDLognormal
        :members:
    .. autoclass:: PSDGatesGaudinSchuhman
        :members:
    .. autoclass:: PSDRosinRammler
        :members:
    .. autoclass:: PSDCustom
        :members:

    Helper functions: Lognormal Distribution
    ----------------------------------------
    .. autofunction:: pdf_lognormal
    .. autofunction:: cdf_lognormal
    .. autofunction:: pdf_lognormal_basis_integral

    Helper functions: Gates Gaudin Schuhman Distribution
    ----------------------------------------------------
    .. autofunction:: pdf_Gates_Gaudin_Schuhman
    .. autofunction:: cdf_Gates_Gaudin_Schuhman
    .. autofunction:: pdf_Gates_Gaudin_Schuhman_basis_integral

    Helper functions: Rosin Rammler Distribution
    --------------------------------------------
    .. autofunction:: pdf_Rosin_Rammler
    .. autofunction:: cdf_Rosin_Rammler
    .. autofunction:: pdf_Rosin_Rammler_basis_integral

    Sieves
    ------
    .. autoclass:: Sieve
    .. autodata:: ASTM_E11_sieves
    .. autodata:: ISO_3310_1_sieves
    .. autodata:: ISO_3310_1_R20
    .. autodata:: ISO_3310_1_R20_3
    .. autodata:: ISO_3310_1_R40_3
    .. autodata:: ISO_3310_1_R10

    Point Spacing
    -------------
    .. autofunction:: psd_spacing
    """

    __all__ = ['ParticleSizeDistribution', 'ParticleSizeDistributionContinuous',
               'PSDLognormal', 'PSDGatesGaudinSchuhman', 'PSDRosinRammler',
               'PSDInterpolated', 'PSDCustom',
               'psd_spacing',

               'pdf_lognormal', 'cdf_lognormal', 'pdf_lognormal_basis_integral',

               'pdf_Gates_Gaudin_Schuhman', 'cdf_Gates_Gaudin_Schuhman',
               'pdf_Gates_Gaudin_Schuhman_basis_integral',

               'pdf_Rosin_Rammler', 'cdf_Rosin_Rammler',
               'pdf_Rosin_Rammler_basis_integral',

               'ASTM_E11_sieves', 'ISO_3310_1_sieves', 'Sieve',
               'ISO_3310_1_R20_3', 'ISO_3310_1_R20', 'ISO_3310_1_R10',
               'ISO_3310_1_R40_3']

    from math import exp, log, log10, pi, sqrt

    from fluids.numerics import brenth, cumsum, diff, epsilon, erf, gamma, gammaincc, linspace, logspace, normalize, quad

    ROOT_TWO_PI = sqrt(2.0*pi)

    NO_MATPLOTLIB_MSG = 'Optional dependency matplotlib is required for plotting'

    class Sieve:
        r'''Class for storing data on sieves. If a property is not available, it is
        set to None.

        Attributes
        ----------
        designation : str
            The standard name of the sieve - its opening's length in units of
            millimeters
        old_designation : str
            The older, imperial-esque name of the sieve; in Numbers, or inches for
            large sieves
        opening : float
            The opening length of the sieve holes, [m]
        opening_inch : float
            The opening length of the sieve holes in the rounded inches as stated
            in common tables (not exactly equal to the `opening`), [inch]
        Y_variation_avg : float
            The allowable average variation in the Y direction of the sieve
            openings, [m]
        X_variation_max : float
            The allowable maximum variation in the X direction of the sieve
            openings, [m]
        max_opening : float
            The maximum allowable opening of the sieve, [m]
        calibration_samples : float
            The number of opening sample inspections required for `calibration`-
            type sieve openings (per 100 ft^2 of sieve material), [1/(ft^2)]
        compliance_sd : float
            The maximum standard deviation of `compliance`-type sieve openings,
            [-]
        inspection_samples : float
            The number of opening sample inspections required for `inspection`-
            type sieve openings (based on an 8-inch sieve), [-]
        inspection_sd : float
            The maximum standard deviation of `inspection`-type sieve openings,
            [-]
        calibration_samples : float
            The number of opening sample inspections required for `calibration`-
            type sieve openings (based on an 8-inch sieve), [-]
        calibration_sd : float
            The maximum standard deviation of `calibration`-type sieve openings,
            [-]
        d_wire : float
            Typical wire diameter of the specified sieve size, [m]
        d_wire_min : float
            Permissible minimum wire diameter of specified sieve size, [m]
        d_wire_max : float
            Permissible maximum wire diameter of specified sieve size, [m]

        '''

        __slots__ = ('designation', 'old_designation', 'opening', 'opening_inch',
                     'Y_variation_avg', 'X_variation_max', 'max_opening',
                     'calibration_samples', 'compliance_sd', 'inspection_samples',
                     'inspection_sd', 'calibration_samples', 'calibration_sd',
                     'd_wire', 'd_wire_min', 'd_wire_max', 'compliance_samples')

    #    def __repr__(self):
    #        s = 'Sieve(%s)'
    #        s2 = ''
    #        for attr, value in self.__dict__.items():
    #            if value is not None:
    #                if type(value) == float:
    #                    value = round(value, 8)
    #                elif type(value) == str:
    #                    value = "'" + value + "'"
    #                s2 += '%s=%s, '%(attr, value)
    #        s2 = s2[0:-2]
    #        return s %(s2)

        def __repr__(self):
            return f'<Sieve, designation {self.designation} mm, opening {self.opening:g} m>'

        def __init__(self, designation, old_designation=None, opening=None,
                     opening_inch=None, Y_variation_avg=None, X_variation_max=None,
                     max_opening=None, compliance_samples=None, compliance_sd=None,
                     inspection_samples=None, inspection_sd=None, calibration_samples=None,
                     calibration_sd=None, d_wire=None, d_wire_min=None, d_wire_max=None):

            self.designation = designation
            self.old_designation = old_designation
            self.opening_inch = opening_inch
            self.opening = opening

            self.Y_variation_avg = Y_variation_avg
            self.X_variation_max = X_variation_max
            self.max_opening = max_opening

            self.compliance_samples = compliance_samples
            self.compliance_sd = compliance_sd

            self.inspection_samples = inspection_samples
            self.inspection_sd = inspection_sd

            self.calibration_samples = calibration_samples
            self.calibration_sd = calibration_sd

            self.d_wire = d_wire
            self.d_wire_min = d_wire_min
            self.d_wire_max = d_wire_max

    ASTM_E11_sieves = {'0.02': Sieve(calibration_samples=300.0, d_wire_min=2e-08, d_wire=2e-08, inspection_sd=4.51, calibration_sd=4.75, old_designation='No. 635', opening=2e-05, compliance_samples=1000.0, opening_inch=8e-07, inspection_samples=100.0, designation='0.02', d_wire_max=2e-08, max_opening=0.035, X_variation_max=1.5e-05, Y_variation_avg=2.3e-06, compliance_sd=5.33),
     '0.025': Sieve(calibration_samples=300.0, d_wire_min=2e-08, d_wire=3e-08, inspection_sd=4.82, calibration_sd=5.06, old_designation='No. 500', opening=2.5e-05, compliance_samples=1000.0, opening_inch=1e-06, inspection_samples=100.0, designation='0.025', d_wire_max=3e-08, max_opening=0.041, X_variation_max=1.6e-05, Y_variation_avg=2.5e-06, compliance_sd=5.71),
     '0.032': Sieve(calibration_samples=300.0, d_wire_min=2e-08, d_wire=3e-08, inspection_sd=5.42, calibration_sd=5.71, old_designation='No. 450', opening=3.2e-05, compliance_samples=1000.0, opening_inch=1.2e-06, inspection_samples=100.0, designation='0.032', d_wire_max=3e-08, max_opening=0.05, X_variation_max=1.8e-05, Y_variation_avg=2.7e-06, compliance_sd=6.42),
     '0.038': Sieve(calibration_samples=300.0, d_wire_min=2e-08, d_wire=3e-08, inspection_sd=5.99, calibration_sd=6.31, old_designation='No. 400', opening=3.8e-05, compliance_samples=1000.0, opening_inch=1.5e-06, inspection_samples=100.0, designation='0.038', d_wire_max=3e-08, max_opening=0.058, X_variation_max=2e-05, Y_variation_avg=2.9e-06, compliance_sd=7.09),
     '0.045': Sieve(calibration_samples=250.0, d_wire_min=3e-08, d_wire=3e-08, inspection_sd=6.56, calibration_sd=6.84, old_designation='No. 325', opening=4.5e-05, compliance_samples=1000.0, opening_inch=1.7e-06, inspection_samples=100.0, designation='0.045', d_wire_max=4e-08, max_opening=0.067, X_variation_max=2.2e-05, Y_variation_avg=3.1e-06, compliance_sd=7.76),
     '0.053': Sieve(calibration_samples=250.0, d_wire_min=3e-08, d_wire=4e-08, inspection_sd=7.13, calibration_sd=7.44, old_designation='No. 270', opening=5.3e-05, compliance_samples=1000.0, opening_inch=2.1e-06, inspection_samples=100.0, designation='0.053', d_wire_max=4e-08, max_opening=0.077, X_variation_max=2.4e-05, Y_variation_avg=3.4e-06, compliance_sd=8.44),
     '0.063': Sieve(calibration_samples=250.0, d_wire_min=4e-08, d_wire=5e-08, inspection_sd=7.76, calibration_sd=8.09, old_designation='No. 230', opening=6.3e-05, compliance_samples=1000.0, opening_inch=2.5e-06, inspection_samples=100.0, designation='0.063', d_wire_max=5e-08, max_opening=0.089, X_variation_max=2.6e-05, Y_variation_avg=3.7e-06, compliance_sd=9.18),
     '0.075': Sieve(calibration_samples=250.0, d_wire_min=4e-08, d_wire=5e-08, inspection_sd=8.64, calibration_sd=9.02, old_designation='No. 200', opening=7.5e-05, compliance_samples=1000.0, opening_inch=2.9e-06, inspection_samples=100.0, designation='0.075', d_wire_max=6e-08, max_opening=0.104, X_variation_max=2.9e-05, Y_variation_avg=4.1e-06, compliance_sd=10.23),
     '0.09': Sieve(calibration_samples=200.0, d_wire_min=5e-08, d_wire=6e-08, inspection_sd=9.53, calibration_sd=9.8, old_designation='No. 170', opening=9e-05, compliance_samples=1000.0, opening_inch=3.5e-06, inspection_samples=100.0, designation='0.09', d_wire_max=7e-08, max_opening=0.122, X_variation_max=3.2e-05, Y_variation_avg=4.6e-06, compliance_sd=11.27),
     '0.106': Sieve(calibration_samples=200.0, d_wire_min=6e-08, d_wire=7e-08, inspection_sd=10.47, calibration_sd=10.77, old_designation='No. 140', opening=0.000106, compliance_samples=1000.0, opening_inch=4.1e-06, inspection_samples=100.0, designation='0.106', d_wire_max=8e-08, max_opening=0.141, X_variation_max=3.5e-05, Y_variation_avg=5.2e-06, compliance_sd=12.39),
     '0.125': Sieve(calibration_samples=200.0, d_wire_min=8e-08, d_wire=9e-08, inspection_sd=11.41, calibration_sd=11.74, old_designation='No. 120', opening=0.000125, compliance_samples=1000.0, opening_inch=4.9e-06, inspection_samples=100.0, designation='0.125', d_wire_max=1e-07, max_opening=0.163, X_variation_max=3.8e-05, Y_variation_avg=5.8e-06, compliance_sd=13.51),
     '0.15': Sieve(calibration_samples=200.0, d_wire_min=9e-08, d_wire=1e-07, inspection_sd=12.93, calibration_sd=13.3, old_designation='No. 100', opening=0.00015, compliance_samples=1000.0, opening_inch=5.9e-06, inspection_samples=100.0, designation='0.15', d_wire_max=1.2e-07, max_opening=0.193, X_variation_max=4.3e-05, Y_variation_avg=6.6e-06, compliance_sd=15.3),
     '0.18': Sieve(calibration_samples=200.0, d_wire_min=1.1e-07, d_wire=1.2e-07, inspection_sd=14.24, calibration_sd=14.65, old_designation='No. 80', opening=0.00018, compliance_samples=1000.0, opening_inch=7e-06, inspection_samples=100.0, designation='0.18', d_wire_max=1.5e-07, max_opening=0.227, X_variation_max=4.7e-05, Y_variation_avg=7.6e-06, compliance_sd=16.85),
     '0.212': Sieve(calibration_samples=160.0, d_wire_min=1.2e-07, d_wire=1.4e-07, inspection_sd=15.59, calibration_sd=16.08, old_designation='No. 70', opening=0.000212, compliance_samples=800.0, opening_inch=8.3e-06, inspection_samples=80.0, designation='0.212', d_wire_max=1.7e-07, max_opening=0.264, X_variation_max=5.2e-05, Y_variation_avg=8.7e-06, compliance_sd=18.79),
     '0.25': Sieve(calibration_samples=160.0, d_wire_min=1.3e-07, d_wire=1.6e-07, inspection_sd=17.44, calibration_sd=17.99, old_designation='No. 60', opening=0.00025, compliance_samples=800.0, opening_inch=9.8e-06, inspection_samples=80.0, designation='0.25', d_wire_max=1.9e-07, max_opening=0.308, X_variation_max=5.8e-05, Y_variation_avg=9.9e-06, compliance_sd=21.02),
     '0.3': Sieve(calibration_samples=160.0, d_wire_min=1.7e-07, d_wire=2e-07, inspection_sd=19.66, calibration_sd=20.29, old_designation='No. 50', opening=0.0003, compliance_samples=800.0, opening_inch=1.17e-05, inspection_samples=80.0, designation='0.3', d_wire_max=2.3e-07, max_opening=0.365, X_variation_max=6.5e-05, Y_variation_avg=1.15e-05, compliance_sd=23.7),
     '0.355': Sieve(calibration_samples=160.0, d_wire_min=1.9e-07, d_wire=2.2e-07, inspection_sd=21.95, calibration_sd=22.64, old_designation='No. 45', opening=0.000355, compliance_samples=800.0, opening_inch=1.39e-05, inspection_samples=80.0, designation='0.355', d_wire_max=2.6e-07, max_opening=0.427, X_variation_max=7.2e-05, Y_variation_avg=1.33e-05, compliance_sd=26.45),
     '0.425': Sieve(calibration_samples=120.0, d_wire_min=2.4e-07, d_wire=2.8e-07, inspection_sd=24.2, calibration_sd=25.08, old_designation='No. 40', opening=0.000425, compliance_samples=600.0, opening_inch=1.65e-05, inspection_samples=60.0, designation='0.425', d_wire_max=3.2e-07, max_opening=0.506, X_variation_max=8.1e-05, Y_variation_avg=1.55e-05, compliance_sd=29.95),
     '0.5': Sieve(calibration_samples=120.0, d_wire_min=2.7e-07, d_wire=3.2e-07, inspection_sd=26.85, calibration_sd=27.82, old_designation='No. 35', opening=0.0005, compliance_samples=600.0, opening_inch=1.97e-05, inspection_samples=60.0, designation='0.5', d_wire_max=3.6e-07, max_opening=0.589, X_variation_max=8.9e-05, Y_variation_avg=1.8e-05, compliance_sd=33.23),
     '0.6': Sieve(calibration_samples=100.0, d_wire_min=3.4e-07, d_wire=4e-07, inspection_sd=30.14, calibration_sd=31.32, old_designation='No. 30', opening=0.0006, compliance_samples=500.0, opening_inch=2.34e-05, inspection_samples=50.0, designation='0.6', d_wire_max=4.6e-07, max_opening=0.701, X_variation_max=0.000101, Y_variation_avg=2.12e-05, compliance_sd=38.0),
     '0.71': Sieve(calibration_samples=100.0, d_wire_min=3.8e-07, d_wire=4.5e-07, inspection_sd=33.82, calibration_sd=35.14, old_designation='No. 25', opening=0.00071, compliance_samples=500.0, opening_inch=2.78e-05, inspection_samples=50.0, designation='0.71', d_wire_max=5.2e-07, max_opening=0.822, X_variation_max=0.000112, Y_variation_avg=2.47e-05, compliance_sd=42.63),
     '0.85': Sieve(calibration_samples=80.0, d_wire_min=4.3e-07, d_wire=5e-07, inspection_sd=37.73, calibration_sd=39.36, old_designation='No. 20', opening=0.00085, compliance_samples=400.0, opening_inch=3.31e-05, inspection_samples=40.0, designation='0.85', d_wire_max=5.8e-07, max_opening=0.977, X_variation_max=0.000127, Y_variation_avg=2.91e-05, compliance_sd=48.76),
     '1': Sieve(calibration_samples=80.0, d_wire_min=0.00048, d_wire=0.00056, inspection_sd=0.042, calibration_sd=0.044, old_designation='No. 18', opening=0.001, compliance_samples=400.0, opening_inch=3.94e-05, inspection_samples=40.0, designation='1', d_wire_max=0.00064, max_opening=1.14, X_variation_max=0.00014, Y_variation_avg=3.4e-05, compliance_sd=0.055),
     '1.18': Sieve(calibration_samples=80.0, d_wire_min=0.00054, d_wire=0.00063, inspection_sd=0.049, calibration_sd=0.051, old_designation='No. 16', opening=0.00118, compliance_samples=400.0, opening_inch=4.69e-05, inspection_samples=40.0, designation='1.18', d_wire_max=0.00072, max_opening=1.34, X_variation_max=0.00016, Y_variation_avg=4e-05, compliance_sd=0.063),
     '1.4': Sieve(calibration_samples=80.0, d_wire_min=0.0006, d_wire=0.00071, inspection_sd=0.055, calibration_sd=0.057, old_designation='No. 14', opening=0.0014, compliance_samples=400.0, opening_inch=5.55e-05, inspection_samples=40.0, designation='1.4', d_wire_max=0.00082, max_opening=1.58, X_variation_max=0.00018, Y_variation_avg=4.6e-05, compliance_sd=0.071),
     '1.7': Sieve(calibration_samples=50.0, d_wire_min=0.00068, d_wire=0.0008, inspection_sd=0.059, calibration_sd=0.062, old_designation='No. 12', opening=0.0017, compliance_samples=250.0, opening_inch=6.61e-05, inspection_samples=25.0, designation='1.7', d_wire_max=0.00092, max_opening=1.9, X_variation_max=0.0002, Y_variation_avg=5.6e-05, compliance_sd=0.081),
     '100': Sieve(d_wire_min=0.0054, d_wire=0.0063, old_designation='4 in.', opening=0.1, compliance_samples=20.0, opening_inch=0.004, designation='100', d_wire_max=0.0072, max_opening=103.82, X_variation_max=0.00382, Y_variation_avg=0.00294),
     '106': Sieve(d_wire_min=0.0054, d_wire=0.0063, old_designation='4.24 in.', opening=0.106, compliance_samples=20.0, opening_inch=0.00424, designation='106', d_wire_max=0.0072, max_opening=109.99, X_variation_max=0.00399, Y_variation_avg=0.00312),
     '11.2': Sieve(calibration_samples=30.0, d_wire_min=0.0021, d_wire=0.0025, inspection_sd=0.256, calibration_sd=0.274, old_designation='7/16 in.', opening=0.0112, compliance_samples=150.0, opening_inch=0.000438, inspection_samples=15.0, designation='11.2', d_wire_max=0.0029, max_opening=11.97, X_variation_max=0.00077, Y_variation_avg=0.000346, compliance_sd=0.382),
     '12.5': Sieve(calibration_samples=30.0, d_wire_min=0.0021, d_wire=0.0025, inspection_sd=0.283, calibration_sd=0.302, old_designation='1/2 in.', opening=0.0125, compliance_samples=150.0, opening_inch=0.0005, inspection_samples=15.0, designation='12.5', d_wire_max=0.0029, max_opening=13.33, X_variation_max=0.00083, Y_variation_avg=0.000385, compliance_sd=0.421),
     '125': Sieve(d_wire_min=0.0068, d_wire=0.008, old_designation='5 in.', opening=0.125, compliance_samples=20.0, opening_inch=0.005, designation='125', d_wire_max=0.0092, max_opening=129.51, X_variation_max=0.00451, Y_variation_avg=0.00366),
     '13.2': Sieve(calibration_samples=30.0, d_wire_min=0.0024, d_wire=0.0028, inspection_sd=0.296, calibration_sd=0.316, old_designation='0.530 in.', opening=0.0132, compliance_samples=150.0, opening_inch=0.00053, inspection_samples=15.0, designation='13.2', d_wire_max=0.0032, max_opening=14.06, X_variation_max=0.00086, Y_variation_avg=0.000406, compliance_sd=0.441),
     '16': Sieve(calibration_samples=30.0, d_wire_min=0.0027, d_wire=0.00315, inspection_sd=0.354, calibration_sd=0.378, old_designation='5/8 in.', opening=0.016, compliance_samples=150.0, opening_inch=0.000625, inspection_samples=15.0, designation='16', d_wire_max=0.0036, max_opening=16.99, X_variation_max=0.00099, Y_variation_avg=0.00049, compliance_sd=0.527),
     '19': Sieve(calibration_samples=30.0, d_wire_min=0.0027, d_wire=0.00315, inspection_sd=0.418, calibration_sd=0.446, old_designation='3/4 in.', opening=0.019, compliance_samples=150.0, opening_inch=0.00075, inspection_samples=15.0, designation='19', d_wire_max=0.0035, max_opening=20.13, X_variation_max=0.00113, Y_variation_avg=0.000579, compliance_sd=0.622),
     '2': Sieve(calibration_samples=50.0, d_wire_min=0.00077, d_wire=0.0009, inspection_sd=0.068, calibration_sd=0.072, old_designation='No. 10', opening=0.002, compliance_samples=250.0, opening_inch=7.87e-05, inspection_samples=25.0, designation='2', d_wire_max=0.00104, max_opening=2.23, X_variation_max=0.00023, Y_variation_avg=6.5e-05, compliance_sd=0.094),
     '2.36': Sieve(calibration_samples=40.0, d_wire_min=0.00085, d_wire=0.001, inspection_sd=0.073, calibration_sd=0.077, old_designation='No. 8', opening=0.00236, compliance_samples=200.0, opening_inch=9.37e-05, inspection_samples=20.0, designation='2.36', d_wire_max=0.00115, max_opening=2.61, X_variation_max=0.00025, Y_variation_avg=7.6e-05, compliance_sd=0.104),
     '2.8': Sieve(calibration_samples=40.0, d_wire_min=0.00095, d_wire=0.00112, inspection_sd=0.085, calibration_sd=0.09, old_designation='No. 7', opening=0.0028, compliance_samples=200.0, opening_inch=0.00011, inspection_samples=20.0, designation='2.8', d_wire_max=0.0013, max_opening=3.09, X_variation_max=0.00029, Y_variation_avg=9e-05, compliance_sd=0.121),
     '22.4': Sieve(d_wire_min=0.003, d_wire=0.00355, inspection_sd=0.493, old_designation='7/8 in.', opening=0.0224, compliance_samples=150.0, opening_inch=0.000875, inspection_samples=15.0, designation='22.4', d_wire_max=0.0041, max_opening=23.67, X_variation_max=0.00127, Y_variation_avg=0.000681, compliance_sd=0.734),
     '25': Sieve(d_wire_min=0.003, d_wire=0.00355, inspection_sd=0.553, old_designation='1.00 in.', opening=0.025, compliance_samples=20.0, opening_inch=0.001, inspection_samples=15.0, designation='25', d_wire_max=0.0041, max_opening=26.38, X_variation_max=0.00138, Y_variation_avg=0.000758, compliance_sd=0.823),
     '26.5': Sieve(d_wire_min=0.003, d_wire=0.00355, inspection_sd=0.584, old_designation='1.06 in.', opening=0.0265, compliance_samples=20.0, opening_inch=0.00106, inspection_samples=15.0, designation='26.5', d_wire_max=0.0041, max_opening=27.94, X_variation_max=0.00144, Y_variation_avg=0.000802, compliance_sd=0.869),
     '3.35': Sieve(calibration_samples=40.0, d_wire_min=0.00106, d_wire=0.00125, inspection_sd=0.097, calibration_sd=0.103, old_designation='No. 6', opening=0.00335, compliance_samples=200.0, opening_inch=0.000132, inspection_samples=20.0, designation='3.35', d_wire_max=0.0015, max_opening=3.67, X_variation_max=0.00032, Y_variation_avg=0.000107, compliance_sd=0.138),
     '31.5': Sieve(d_wire_min=0.0034, d_wire=0.004, old_designation='1 1/4 in.', opening=0.0315, compliance_samples=20.0, opening_inch=0.00125, designation='31.5', d_wire_max=0.0046, max_opening=33.13, X_variation_max=0.00163, Y_variation_avg=0.00095, compliance_sd=1.066),
     '37.5': Sieve(d_wire_min=0.0038, d_wire=0.0045, old_designation='1 1/2 in.', opening=0.0375, compliance_samples=20.0, opening_inch=0.0015, designation='37.5', d_wire_max=0.0052, max_opening=39.35, X_variation_max=0.00185, Y_variation_avg=0.00113, compliance_sd=1.374),
     '4': Sieve(calibration_samples=30.0, d_wire_min=0.0012, d_wire=0.0014, inspection_sd=0.108, calibration_sd=0.115, old_designation='No. 5', opening=0.004, compliance_samples=150.0, opening_inch=0.000157, inspection_samples=15.0, designation='4', d_wire_max=0.0017, max_opening=4.37, X_variation_max=0.00037, Y_variation_avg=0.000127, compliance_sd=0.161),
     '4.75': Sieve(calibration_samples=30.0, d_wire_min=0.0013, d_wire=0.0016, inspection_sd=0.123, calibration_sd=0.131, old_designation='No. 4', opening=0.00475, compliance_samples=150.0, opening_inch=0.000187, inspection_samples=15.0, designation='4.75', d_wire_max=0.0019, max_opening=5.16, X_variation_max=0.00041, Y_variation_avg=0.00015, compliance_sd=0.182),
     '45': Sieve(d_wire_min=0.0038, d_wire=0.0045, old_designation='1 3/4 in.', opening=0.045, compliance_samples=20.0, opening_inch=0.00175, designation='45', d_wire_max=0.0052, max_opening=47.12, X_variation_max=0.00212, Y_variation_avg=0.00135),
     '5.6': Sieve(calibration_samples=30.0, d_wire_min=0.0013, d_wire=0.0016, inspection_sd=0.142, calibration_sd=0.151, old_designation='No. 3 1/2', opening=0.0056, compliance_samples=150.0, opening_inch=0.000223, inspection_samples=15.0, designation='5.6', d_wire_max=0.0019, max_opening=6.07, X_variation_max=0.00047, Y_variation_avg=0.000176, compliance_sd=0.211),
     '50': Sieve(d_wire_min=0.0043, d_wire=0.005, old_designation='2 in.', opening=0.05, compliance_samples=20.0, opening_inch=0.002, designation='50', d_wire_max=0.0058, max_opening=52.29, X_variation_max=0.00229, Y_variation_avg=0.00149),
     '53': Sieve(d_wire_min=0.0043, d_wire=0.005, old_designation='2.12 in.', opening=0.053, compliance_samples=20.0, opening_inch=0.00212, designation='53', d_wire_max=0.0058, max_opening=55.39, X_variation_max=0.00239, Y_variation_avg=0.00158),
     '6.3': Sieve(calibration_samples=30.0, d_wire_min=0.0015, d_wire=0.0018, inspection_sd=0.157, calibration_sd=0.167, old_designation='1/4 in.', opening=0.0063, compliance_samples=150.0, opening_inch=0.00025, inspection_samples=15.0, designation='6.3', d_wire_max=0.0021, max_opening=6.81, X_variation_max=0.00051, Y_variation_avg=0.000197, compliance_sd=0.233),
     '6.7': Sieve(calibration_samples=30.0, d_wire_min=0.0015, d_wire=0.0018, inspection_sd=0.164, calibration_sd=0.175, old_designation='0.265 in.', opening=0.0067, compliance_samples=150.0, opening_inch=0.000265, inspection_samples=15.0, designation='6.7', d_wire_max=0.0021, max_opening=7.23, X_variation_max=0.00053, Y_variation_avg=0.00021, compliance_sd=0.245),
     '63': Sieve(d_wire_min=0.0048, d_wire=0.0056, old_designation='2 1/2 in.', opening=0.063, compliance_samples=20.0, opening_inch=0.0025, designation='63', d_wire_max=0.0064, max_opening=65.71, X_variation_max=0.00271, Y_variation_avg=0.00187),
     '75': Sieve(d_wire_min=0.0054, d_wire=0.0063, old_designation='3 in.', opening=0.075, compliance_samples=20.0, opening_inch=0.003, designation='75', d_wire_max=0.0072, max_opening=78.09, X_variation_max=0.00309, Y_variation_avg=0.00222),
     '8': Sieve(calibration_samples=30.0, d_wire_min=0.0017, d_wire=0.002, inspection_sd=0.191, calibration_sd=0.204, old_designation='5/16 in.', opening=0.008, compliance_samples=150.0, opening_inch=0.000312, inspection_samples=15.0, designation='8', d_wire_max=0.0023, max_opening=8.6, X_variation_max=0.0006, Y_variation_avg=0.000249, compliance_sd=0.284),
     '9.5': Sieve(calibration_samples=30.0, d_wire_min=0.0019, d_wire=0.00224, inspection_sd=0.222, calibration_sd=0.237, old_designation='3/8 in.', opening=0.0095, compliance_samples=150.0, opening_inch=0.000375, inspection_samples=15.0, designation='9.5', d_wire_max=0.0026, max_opening=10.18, X_variation_max=0.00068, Y_variation_avg=0.000295, compliance_sd=0.33),
     '90': Sieve(d_wire_min=0.0054, d_wire=0.0063, old_designation='3 1/2 in.', opening=0.09, compliance_samples=20.0, opening_inch=0.0035, designation='90', d_wire_max=0.0072, max_opening=93.53, X_variation_max=0.00353, Y_variation_avg=0.00265)
     }
    """Dictionary containing ASTM E-11 sieve series :py:func:`Sieve` objects, indexed by
    their size in mm as a string.

    References
    ----------
    .. [1] ASTM E11 - 17 - Standard Specification for Woven Wire Test Sieve
       Cloth and Test Sieves.
    """

    ASTM_E11_sieve_designations = ['125', '106', '100', '90', '75', '63', '53',
                                   '50', '45', '37.5', '31.5', '26.5', '25',
                                   '22.4', '19', '16', '13.2', '12.5', '11.2',
                                   '9.5', '8', '6.7', '6.3', '5.6', '4.75', '4',
                                   '3.35', '2.8', '2.36', '2', '1.7', '1.4',
                                   '1.18', '1', '0.85', '0.71', '0.6', '0.5',
                                   '0.425', '0.355', '0.3', '0.25', '0.212',
                                   '0.18', '0.15', '0.125', '0.106', '0.09',
                                   '0.075', '0.063', '0.053', '0.045', '0.038',
                                   '0.032', '0.025', '0.02']

    ASTM_E11_sieve_list = [ASTM_E11_sieves[i] for i in ASTM_E11_sieve_designations]

    ISO_3310_1_sieves = {
     '0.02': Sieve(designation='0.02', d_wire_max=2.3e-05, compliance_sd=4.7e-06, X_variation_max=1.3e-05, d_wire=2e-05, Y_variation_avg=2.1e-06, d_wire_min=2.3e-05, opening=2e-05),
     '0.025': Sieve(designation='0.025', d_wire_max=2.9e-05, compliance_sd=5.2e-06, X_variation_max=1.5e-05, d_wire=2.5e-05, Y_variation_avg=2.2e-06, d_wire_min=2.9e-05, opening=2.5e-05),
     '0.032': Sieve(designation='0.032', d_wire_max=3.3e-05, compliance_sd=5.9e-06, X_variation_max=1.7e-05, d_wire=2.8e-05, Y_variation_avg=2.4e-06, d_wire_min=3.3e-05, opening=3.2e-05),
     '0.036': Sieve(designation='0.036', d_wire_max=3.5e-05, compliance_sd=6.3e-06, X_variation_max=1.8e-05, d_wire=3e-05, Y_variation_avg=2.6e-06, d_wire_min=3.5e-05, opening=3.6e-05),
     '0.038': Sieve(designation='0.038', d_wire_max=3.5e-05, compliance_sd=6.4e-06, X_variation_max=1.8e-05, d_wire=3e-05, Y_variation_avg=2.6e-06, d_wire_min=3.5e-05, opening=3.8e-05),
     '0.04': Sieve(designation='0.04', d_wire_max=3.7e-05, compliance_sd=6.5e-06, X_variation_max=1.9e-05, d_wire=3.2e-05, Y_variation_avg=2.7e-06, d_wire_min=3.7e-05, opening=4e-05),
     '0.045': Sieve(designation='0.045', d_wire_max=3.7e-05, compliance_sd=6.9e-06, X_variation_max=2e-05, d_wire=3.2e-05, Y_variation_avg=2.8e-06, d_wire_min=3.7e-05, opening=4.5e-05),
     '0.05': Sieve(designation='0.05', d_wire_max=4.1e-05, compliance_sd=7.3e-06, X_variation_max=2.1e-05, d_wire=3.6e-05, Y_variation_avg=3e-06, d_wire_min=4.1e-05, opening=5e-05),
     '0.053': Sieve(designation='0.053', d_wire_max=4.1e-05, compliance_sd=7.6e-06, X_variation_max=2.1e-05, d_wire=3.6e-05, Y_variation_avg=3.1e-06, d_wire_min=4.1e-05, opening=5.3e-05),
     '0.056': Sieve(designation='0.056', d_wire_max=4.6e-05, compliance_sd=7.8e-06, X_variation_max=2.2e-05, d_wire=4e-05, Y_variation_avg=3.2e-06, d_wire_min=4.6e-05, opening=5.6e-05),
     '0.063': Sieve(designation='0.063', d_wire_max=5.2e-05, compliance_sd=8.3e-06, X_variation_max=2.4e-05, d_wire=4.5e-05, Y_variation_avg=3.4e-06, d_wire_min=5.2e-05, opening=6.3e-05),
     '0.071': Sieve(designation='0.071', d_wire_max=5.8e-05, compliance_sd=8.9e-06, X_variation_max=2.5e-05, d_wire=5e-05, Y_variation_avg=3.6e-06, d_wire_min=4.3e-05, opening=7.1e-05),
     '0.075': Sieve(designation='0.075', d_wire_max=5.8e-05, compliance_sd=9.1e-06, X_variation_max=2.6e-05, d_wire=5e-05, Y_variation_avg=3.7e-06, d_wire_min=4.3e-05, opening=7.5e-05),
     '0.08': Sieve(designation='0.08', d_wire_max=6.4e-05, compliance_sd=9.4e-06, X_variation_max=2.7e-05, d_wire=5.6e-05, Y_variation_avg=3.9e-06, d_wire_min=4.8e-05, opening=8e-05),
     '0.09': Sieve(designation='0.09', d_wire_max=7.2e-05, compliance_sd=1.01e-05, X_variation_max=2.9e-05, d_wire=6.3e-05, Y_variation_avg=4.2e-06, d_wire_min=5.4e-05, opening=9e-05),
     '0.1': Sieve(designation='0.1', d_wire_max=8.2e-05, compliance_sd=1.08e-05, X_variation_max=3e-05, d_wire=7.1e-05, Y_variation_avg=4.5e-06, d_wire_min=6e-05, opening=0.0001),
     '0.106': Sieve(designation='0.106', d_wire_max=8.2e-05, compliance_sd=1.11e-05, X_variation_max=3.1e-05, d_wire=7.1e-05, Y_variation_avg=4.7e-06, d_wire_min=6e-05, opening=0.000106),
     '0.112': Sieve(designation='0.112', d_wire_max=9.2e-05, compliance_sd=1.15e-05, X_variation_max=3.2e-05, d_wire=8e-05, Y_variation_avg=4.8e-06, d_wire_min=6.8e-05, opening=0.000112),
     '0.125': Sieve(designation='0.125', d_wire_max=0.000104, compliance_sd=1.22e-05, X_variation_max=3.4e-05, d_wire=9e-05, Y_variation_avg=5.2e-06, d_wire_min=7.7e-05, opening=0.000125),
     '0.14': Sieve(designation='0.14', d_wire_max=0.000115, compliance_sd=1.31e-05, X_variation_max=3.7e-05, d_wire=0.0001, Y_variation_avg=5.7e-06, d_wire_min=8.5e-05, opening=0.00014),
     '0.15': Sieve(designation='0.15', d_wire_max=0.000115, compliance_sd=1.37e-05, X_variation_max=3.8e-05, d_wire=0.0001, Y_variation_avg=6e-06, d_wire_min=8.5e-05, opening=0.00015),
     '0.16': Sieve(designation='0.16', d_wire_max=0.00013, compliance_sd=1.42e-05, X_variation_max=4e-05, d_wire=0.000112, Y_variation_avg=6.3e-06, d_wire_min=9.5e-05, opening=0.00016),
     '0.18': Sieve(designation='0.18', d_wire_max=0.00015, compliance_sd=1.53e-05, X_variation_max=4.3e-05, d_wire=0.000125, Y_variation_avg=6.8e-06, d_wire_min=0.000106, opening=0.00018),
     '0.2': Sieve(designation='0.2', d_wire_max=0.00017, compliance_sd=1.63e-05, X_variation_max=4.5e-05, d_wire=0.00014, Y_variation_avg=7.4e-06, d_wire_min=0.00012, opening=0.0002),
     '0.212': Sieve(designation='0.212', d_wire_max=0.00017, compliance_sd=1.69e-05, X_variation_max=4.7e-05, d_wire=0.00014, Y_variation_avg=7.8e-06, d_wire_min=0.00012, opening=0.000212),
     '0.224': Sieve(designation='0.224', d_wire_max=0.00019, compliance_sd=1.75e-05, X_variation_max=4.9e-05, d_wire=0.00016, Y_variation_avg=8.1e-06, d_wire_min=0.00013, opening=0.000224),
     '0.25': Sieve(designation='0.25', d_wire_max=0.00019, compliance_sd=1.88e-05, X_variation_max=5.2e-05, d_wire=0.00016, Y_variation_avg=8.9e-06, d_wire_min=0.00013, opening=0.00025),
     '0.28': Sieve(designation='0.28', d_wire_max=0.00021, compliance_sd=2.03e-05, X_variation_max=5.6e-05, d_wire=0.00018, Y_variation_avg=1e-05, d_wire_min=0.00015, opening=0.00028),
     '0.3': Sieve(designation='0.3', d_wire_max=0.00023, compliance_sd=2.12e-05, X_variation_max=5.8e-05, d_wire=0.0002, Y_variation_avg=1e-05, d_wire_min=0.00017, opening=0.0003),
     '0.315': Sieve(designation='0.315', d_wire_max=0.00023, compliance_sd=2.19e-05, X_variation_max=6e-05, d_wire=0.0002, Y_variation_avg=1.1e-05, d_wire_min=0.00017, opening=0.000315),
     '0.355': Sieve(designation='0.355', d_wire_max=0.00026, compliance_sd=2.37e-05, X_variation_max=6.5e-05, d_wire=0.000224, Y_variation_avg=1.2e-05, d_wire_min=0.00019, opening=0.000355),
     '0.4': Sieve(designation='0.4', d_wire_max=0.00029, compliance_sd=2.57e-05, X_variation_max=7e-05, d_wire=0.00025, Y_variation_avg=1.3e-05, d_wire_min=0.00021, opening=0.0004),
     '0.425': Sieve(designation='0.425', d_wire_max=0.00032, compliance_sd=2.68e-05, X_variation_max=7.3e-05, d_wire=0.00028, Y_variation_avg=1.4e-05, d_wire_min=0.00024, opening=0.000425),
     '0.45': Sieve(designation='0.45', d_wire_max=0.00032, compliance_sd=2.79e-05, X_variation_max=7.5e-05, d_wire=0.00028, Y_variation_avg=1.5e-05, d_wire_min=0.00024, opening=0.00045),
     '0.5': Sieve(designation='0.5', d_wire_max=0.00036, compliance_sd=3e-05, X_variation_max=8e-05, d_wire=0.000315, Y_variation_avg=1.6e-05, d_wire_min=0.00027, opening=0.0005),
     '0.56': Sieve(designation='0.56', d_wire_max=0.00041, compliance_sd=3.24e-05, X_variation_max=8.7e-05, d_wire=0.000355, Y_variation_avg=1.8e-05, d_wire_min=0.0003, opening=0.00056),
     '0.6': Sieve(designation='0.6', d_wire_max=0.00046, compliance_sd=3.4e-05, X_variation_max=9.1e-05, d_wire=0.0004, Y_variation_avg=1.9e-05, d_wire_min=0.00034, opening=0.0006),
     '0.63': Sieve(designation='0.63', d_wire_max=0.00046, compliance_sd=3.52e-05, X_variation_max=9.3e-05, d_wire=0.0004, Y_variation_avg=2e-05, d_wire_min=0.00034, opening=0.00063),
     '0.71': Sieve(designation='0.71', d_wire_max=0.00052, compliance_sd=3.84e-05, X_variation_max=0.000101, d_wire=0.00045, Y_variation_avg=2.2e-05, d_wire_min=0.00038, opening=0.00071),
     '0.8': Sieve(designation='0.8', d_wire_max=0.00052, compliance_sd=4.18e-05, X_variation_max=0.000109, d_wire=0.00045, Y_variation_avg=2.5e-05, d_wire_min=0.00038, opening=0.0008),
     '0.85': Sieve(designation='0.85', d_wire_max=0.00058, compliance_sd=4.36e-05, X_variation_max=0.000114, d_wire=0.0005, Y_variation_avg=2.6e-05, d_wire_min=0.00043, opening=0.00085),
     '0.9': Sieve(designation='0.9', d_wire_max=0.00058, compliance_sd=4.55e-05, X_variation_max=0.000118, d_wire=0.0005, Y_variation_avg=2.8e-05, d_wire_min=0.00043, opening=0.0009),
     '1': Sieve(designation='1', d_wire_max=0.00064, compliance_sd=4.9e-05, X_variation_max=0.00013, d_wire=0.00056, Y_variation_avg=3e-05, d_wire_min=0.00048, opening=0.001),
     '1.12': Sieve(designation='1.12', d_wire_max=0.00064, compliance_sd=5.3e-05, X_variation_max=0.00014, d_wire=0.00056, Y_variation_avg=3e-05, d_wire_min=0.00048, opening=0.00112),
     '1.18': Sieve(designation='1.18', d_wire_max=0.00072, compliance_sd=5.6e-05, X_variation_max=0.00014, d_wire=0.00063, Y_variation_avg=4e-05, d_wire_min=0.00054, opening=0.00118),
     '1.25': Sieve(designation='1.25', d_wire_max=0.00072, compliance_sd=5.8e-05, X_variation_max=0.00015, d_wire=0.00063, Y_variation_avg=4e-05, d_wire_min=0.00054, opening=0.00125),
     '1.4': Sieve(designation='1.4', d_wire_max=0.00082, compliance_sd=6.3e-05, X_variation_max=0.00016, d_wire=0.00071, Y_variation_avg=4e-05, d_wire_min=0.0006, opening=0.0014),
     '1.6': Sieve(designation='1.6', d_wire_max=0.00092, compliance_sd=7e-05, X_variation_max=0.00017, d_wire=0.0008, Y_variation_avg=5e-05, d_wire_min=0.00068, opening=0.0016),
     '1.7': Sieve(designation='1.7', d_wire_max=0.00092, compliance_sd=7.3e-05, X_variation_max=0.00018, d_wire=0.0008, Y_variation_avg=5e-05, d_wire_min=0.00068, opening=0.0017),
     '1.8': Sieve(designation='1.8', d_wire_max=0.00092, compliance_sd=7.6e-05, X_variation_max=0.00019, d_wire=0.0008, Y_variation_avg=5e-05, d_wire_min=0.00068, opening=0.0018),
     '10': Sieve(designation='10', d_wire_max=0.0029, compliance_sd=0.000307, X_variation_max=0.00064, d_wire=0.0025, Y_variation_avg=0.00028, d_wire_min=0.0021, opening=0.01),
     '100': Sieve(designation='100', d_wire_max=0.0072, X_variation_max=0.00344, d_wire=0.0063, Y_variation_avg=0.00265, d_wire_min=0.0054, opening=0.1),
     '106': Sieve(designation='106', d_wire_max=0.0072, X_variation_max=0.00359, d_wire=0.0063, Y_variation_avg=0.0028, d_wire_min=0.0054, opening=0.106),
     '11.2': Sieve(designation='11.2', d_wire_max=0.0029, compliance_sd=0.000339, X_variation_max=0.00069, d_wire=0.0025, Y_variation_avg=0.00031, d_wire_min=0.0021, opening=0.0112),
     '112': Sieve(designation='112', d_wire_max=0.0092, X_variation_max=0.00374, d_wire=0.008, Y_variation_avg=0.00296, d_wire_min=0.0068, opening=0.112),
     '12.5': Sieve(designation='12.5', d_wire_max=0.0029, compliance_sd=0.000374, X_variation_max=0.00075, d_wire=0.0025, Y_variation_avg=0.00035, d_wire_min=0.0021, opening=0.0125),
     '125': Sieve(designation='125', d_wire_max=0.0092, X_variation_max=0.00406, d_wire=0.008, Y_variation_avg=0.0033, d_wire_min=0.0068, opening=0.125),
     '13.2': Sieve(designation='13.2', d_wire_max=0.0032, compliance_sd=0.000392, X_variation_max=0.00078, d_wire=0.0028, Y_variation_avg=0.00037, d_wire_min=0.0024, opening=0.0132),
     '14': Sieve(designation='14', d_wire_max=0.0032, compliance_sd=0.000413, X_variation_max=0.00081, d_wire=0.0028, Y_variation_avg=0.00039, d_wire_min=0.0024, opening=0.014),
     '16': Sieve(designation='16', d_wire_max=0.0036, compliance_sd=0.000467, X_variation_max=0.00089, d_wire=0.00315, Y_variation_avg=0.00044, d_wire_min=0.0027, opening=0.016),
     '18': Sieve(designation='18', d_wire_max=0.0036, compliance_sd=0.00052, X_variation_max=0.00097, d_wire=0.00315, Y_variation_avg=0.00049, d_wire_min=0.0027, opening=0.018),
     '19': Sieve(designation='19', d_wire_max=0.0036, compliance_sd=0.000548, X_variation_max=0.00101, d_wire=0.00315, Y_variation_avg=0.00052, d_wire_min=0.0027, opening=0.019),
     '2': Sieve(designation='2', d_wire_max=0.00104, compliance_sd=8.3e-05, X_variation_max=0.0002, d_wire=0.0009, Y_variation_avg=6e-05, d_wire_min=0.00077, opening=0.002),
     '2.24': Sieve(designation='2.24', d_wire_max=0.00104, compliance_sd=9e-05, X_variation_max=0.00022, d_wire=0.0009, Y_variation_avg=7e-05, d_wire_min=0.00077, opening=0.00224),
     '2.36': Sieve(designation='2.36', d_wire_max=0.00115, compliance_sd=9.4e-05, X_variation_max=0.00023, d_wire=0.001, Y_variation_avg=7e-05, d_wire_min=0.00085, opening=0.00236),
     '2.5': Sieve(designation='2.5', d_wire_max=0.00115, compliance_sd=9.8e-05, X_variation_max=0.00024, d_wire=0.001, Y_variation_avg=7e-05, d_wire_min=0.00085, opening=0.0025),
     '2.8': Sieve(designation='2.8', d_wire_max=0.0013, compliance_sd=0.000108, X_variation_max=0.00026, d_wire=0.00112, Y_variation_avg=8e-05, d_wire_min=0.00095, opening=0.0028),
     '20': Sieve(designation='20', d_wire_max=0.0036, compliance_sd=0.000575, X_variation_max=0.00105, d_wire=0.00315, Y_variation_avg=0.00055, d_wire_min=0.0027, opening=0.02),
     '22.4': Sieve(designation='22.4', d_wire_max=0.0041, compliance_sd=0.000641, X_variation_max=0.00114, d_wire=0.00355, Y_variation_avg=0.00061, d_wire_min=0.003, opening=0.0224),
     '25': Sieve(designation='25', d_wire_max=0.0041, compliance_sd=0.000713, X_variation_max=0.00124, d_wire=0.00355, Y_variation_avg=0.00068, d_wire_min=0.003, opening=0.025),
     '26.5': Sieve(designation='26.5', d_wire_max=0.0041, compliance_sd=0.000757, X_variation_max=0.00129, d_wire=0.00355, Y_variation_avg=0.00072, d_wire_min=0.003, opening=0.0265),
     '28': Sieve(designation='28', d_wire_max=0.0041, compliance_sd=0.000801, X_variation_max=0.00135, d_wire=0.00355, Y_variation_avg=0.00076, d_wire_min=0.003, opening=0.028),
     '3.15': Sieve(designation='3.15', d_wire_max=0.0015, compliance_sd=0.000118, X_variation_max=0.00028, d_wire=0.00125, Y_variation_avg=9e-05, d_wire_min=0.00106, opening=0.00315),
     '3.35': Sieve(designation='3.35', d_wire_max=0.0015, compliance_sd=0.000124, X_variation_max=0.00029, d_wire=0.00125, Y_variation_avg=0.0001, d_wire_min=0.00106, opening=0.00335),
     '3.55': Sieve(designation='3.55', d_wire_max=0.0015, compliance_sd=0.00013, X_variation_max=0.0003, d_wire=0.00125, Y_variation_avg=0.0001, d_wire_min=0.00106, opening=0.00355),
     '31.5': Sieve(designation='31.5', d_wire_max=0.0046, compliance_sd=0.000905, X_variation_max=0.00147, d_wire=0.004, Y_variation_avg=0.00085, d_wire_min=0.0034, opening=0.0315),
     '35.5': Sieve(designation='35.5', d_wire_max=0.0046, compliance_sd=0.001, X_variation_max=0.0016, d_wire=0.004, Y_variation_avg=0.00096, d_wire_min=0.0034, opening=0.0355),
     '37.5': Sieve(designation='37.5', d_wire_max=0.0052, compliance_sd=0.001, X_variation_max=0.00167, d_wire=0.0045, Y_variation_avg=0.00101, d_wire_min=0.0038, opening=0.0375),
     '4': Sieve(designation='4', d_wire_max=0.0017, compliance_sd=0.000143, X_variation_max=0.00033, d_wire=0.0014, Y_variation_avg=0.00011, d_wire_min=0.0012, opening=0.004),
     '4.5': Sieve(designation='4.5', d_wire_max=0.0017, compliance_sd=0.000157, X_variation_max=0.00036, d_wire=0.0014, Y_variation_avg=0.00013, d_wire_min=0.0012, opening=0.0045),
     '4.75': Sieve(designation='4.75', d_wire_max=0.0019, compliance_sd=0.000164, X_variation_max=0.00037, d_wire=0.0016, Y_variation_avg=0.00014, d_wire_min=0.0013, opening=0.00475),
     '40': Sieve(designation='40', d_wire_max=0.0052, compliance_sd=0.001, X_variation_max=0.00175, d_wire=0.0045, Y_variation_avg=0.00108, d_wire_min=0.0038, opening=0.04),
     '45': Sieve(designation='45', d_wire_max=0.0052, compliance_sd=0.001, X_variation_max=0.00191, d_wire=0.0045, Y_variation_avg=0.00121, d_wire_min=0.0038, opening=0.045),
     '5': Sieve(designation='5', d_wire_max=0.0019, compliance_sd=0.000171, X_variation_max=0.00039, d_wire=0.0016, Y_variation_avg=0.00014, d_wire_min=0.0013, opening=0.005),
     '5.6': Sieve(designation='5.6', d_wire_max=0.0019, compliance_sd=0.000188, X_variation_max=0.00042, d_wire=0.0016, Y_variation_avg=0.00016, d_wire_min=0.0013, opening=0.0056),
     '50': Sieve(designation='50', d_wire_max=0.0058, X_variation_max=0.00206, d_wire=0.005, Y_variation_avg=0.00134, d_wire_min=0.0043, opening=0.05),
     '53': Sieve(designation='53', d_wire_max=0.0058, X_variation_max=0.00215, d_wire=0.005, Y_variation_avg=0.00142, d_wire_min=0.0043, opening=0.053),
     '56': Sieve(designation='56', d_wire_max=0.0058, X_variation_max=0.00224, d_wire=0.005, Y_variation_avg=0.0015, d_wire_min=0.0043, opening=0.056),
     '6.3': Sieve(designation='6.3', d_wire_max=0.0021, compliance_sd=0.000207, X_variation_max=0.00046, d_wire=0.0018, Y_variation_avg=0.00018, d_wire_min=0.0015, opening=0.0063),
     '6.7': Sieve(designation='6.7', d_wire_max=0.0021, compliance_sd=0.000218, X_variation_max=0.00048, d_wire=0.0018, Y_variation_avg=0.00019, d_wire_min=0.0015, opening=0.0067),
     '63': Sieve(designation='63', d_wire_max=0.0064, X_variation_max=0.00244, d_wire=0.0056, Y_variation_avg=0.00169, d_wire_min=0.0048, opening=0.063),
     '7.1': Sieve(designation='7.1', d_wire_max=0.0021, compliance_sd=0.000229, X_variation_max=0.0005, d_wire=0.0018, Y_variation_avg=0.0002, d_wire_min=0.0015, opening=0.0071),
     '71': Sieve(designation='71', d_wire_max=0.0064, X_variation_max=0.00267, d_wire=0.0056, Y_variation_avg=0.00189, d_wire_min=0.0048, opening=0.071),
     '75': Sieve(designation='75', d_wire_max=0.0072, X_variation_max=0.00278, d_wire=0.0063, Y_variation_avg=0.002, d_wire_min=0.0054, opening=0.075),
     '8': Sieve(designation='8', d_wire_max=0.0023, compliance_sd=0.000254, X_variation_max=0.00054, d_wire=0.002, Y_variation_avg=0.00022, d_wire_min=0.0017, opening=0.008),
     '80': Sieve(designation='80', d_wire_max=0.0072, X_variation_max=0.00291, d_wire=0.0063, Y_variation_avg=0.00213, d_wire_min=0.0054, opening=0.08),
     '9': Sieve(designation='9', d_wire_max=0.0026, compliance_sd=0.000281, X_variation_max=0.00059, d_wire=0.00224, Y_variation_avg=0.00025, d_wire_min=0.0019, opening=0.009),
     '9.5': Sieve(designation='9.5', d_wire_max=0.0026, compliance_sd=0.000294, X_variation_max=0.00061, d_wire=0.00224, Y_variation_avg=0.00027, d_wire_min=0.0019, opening=0.0095),
     '90': Sieve(designation='90', d_wire_max=0.0072, X_variation_max=0.00318, d_wire=0.0063, Y_variation_avg=0.00239, d_wire_min=0.0054, opening=0.09)
    }
    """Dictionary containing all of the individual :py:func:`Sieve` objects, on the
    ISO 3310-1:2016 series, indexed by their size in mm as a string.

    References
    ----------
    .. [1] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
       -- Part 1: Test Sieves of Metal Wire Cloth.
    """

    ISO_3310_1_sieve_designations  = ['125', '112', '106', '100', '90', '80', '75',
                                      '71', '63', '56', '53', '50', '45', '40',
                                      '37.5', '35.5', '31.5', '28', '26.5', '25',
                                      '22.4', '20', '19', '18', '16', '14', '13.2',
                                      '12.5', '11.2', '10', '9.5', '9', '8', '7.1',
                                      '6.7', '6.3', '5.6', '5', '4.75', '4.5', '4',
                                      '3.55', '3.35', '3.15', '2.8', '2.5', '2.36',
                                      '2.24', '2', '1.8', '1.7', '1.6', '1.4',
                                      '1.25', '1.18', '1.12', '1', '0.9', '0.85',
                                      '0.8', '0.71', '0.63', '0.6', '0.56', '0.5',
                                      '0.45', '0.425', '0.4', '0.355', '0.315',
                                      '0.3', '0.28', '0.25', '0.224', '0.212',
                                      '0.2', '0.18', '0.16', '0.15', '0.14',
                                      '0.125', '0.112', '0.106', '0.1', '0.09',
                                      '0.08', '0.075', '0.071', '0.063', '0.056',
                                      '0.053', '0.05', '0.045', '0.04', '0.038',
                                      '0.036', '0.032', '0.025', '0.02']

    ISO_3310_1_sieve_list = [ISO_3310_1_sieves[i] for i in ISO_3310_1_sieve_designations]

    ISO_3310_1_R20_3 = [ISO_3310_1_sieves[i] for i in ('125', '90', '63', '45', '31.5', '22.4', '16', '11.2', '8', '5.6', '4', '2.8', '2', '1.4', '1', '0.71', '0.5', '0.355', '0.25', '0.18', '0.125', '0.09', '0.063', '0.045')]
    """List containing all of the individual :py:func:`Sieve` objects, on the
    ISO 3310-1:2016 R20/3 series only, ordered from largest  openings to smallest.

    References
    ----------
    .. [1] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
       -- Part 1: Test Sieves of Metal Wire Cloth.
    """

    ISO_3310_1_R20 = ['125', '112', '100', '90', '80', '71', '63', '56', '50', '45', '40', '35.5', '31.5', '28', '25', '22.4', '20', '18', '16', '14', '12.5', '11.2', '10', '9', '8', '7.1', '6.3', '5.6', '5', '4.5', '4', '3.55', '3.15', '2.8', '2.5', '2.24', '2', '1.8', '1.6', '1.4', '1.25', '1.12', '1', '0.9', '0.8', '0.71', '0.63', '0.56', '0.5', '0.45', '0.4', '0.355', '0.315', '0.28', '0.25', '0.224', '0.2', '0.18', '0.16', '0.14', '0.125', '0.112', '0.1', '0.09', '0.08', '0.071', '0.063', '0.056', '0.05', '0.045', '0.04', '0.036']
    ISO_3310_1_R20 = [ISO_3310_1_sieves[i] for i in ISO_3310_1_R20]
    """List containing all of the individual :py:func:`Sieve` objects, on the
    ISO 3310-1:2016 R20 series only, ordered from largest  openings to smallest.

    References
    ----------
    .. [1] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
       -- Part 1: Test Sieves of Metal Wire Cloth.
    """

    ISO_3310_1_R40_3 = ['125', '106', '90', '75', '63', '53', '45', '37.5', '31.5', '26.5', '22.4', '19', '16', '13.2', '11.2', '9.5', '8', '6.7', '5.6', '4.75', '4', '3.35', '2.8', '2.36', '2', '1.7', '1.4', '1.18', '1', '0.85', '0.71', '0.6', '0.5', '0.425', '0.355', '0.3', '0.25', '0.212', '0.18', '0.15', '0.125', '0.106', '0.09', '0.075', '0.063', '0.053', '0.045', '0.038']
    ISO_3310_1_R40_3 = [ISO_3310_1_sieves[i] for i in ISO_3310_1_R40_3]
    """List containing all of the individual :py:func:`Sieve` objects, on the
    ISO 3310-1:2016 R40/3 series only, ordered from largest  openings to smallest.

    References
    ----------
    .. [1] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
       -- Part 1: Test Sieves of Metal Wire Cloth.
    """

    ISO_3310_1_R10 = ['0.036', '0.032', '0.025', '0.02']
    ISO_3310_1_R10 = [ISO_3310_1_sieves[i] for i in ISO_3310_1_R10]
    """List containing all of the individual :py:func:`Sieve` objects, on the
    ISO 3310-1:2016 R10 series only, ordered from largest  openings to smallest.

    References
    ----------
    .. [1] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
       -- Part 1: Test Sieves of Metal Wire Cloth.
    """

    sieve_spacing_options = {'ISO 3310-1': ISO_3310_1_sieve_list,
                             'ISO 3310-1 R20': ISO_3310_1_R20,
                             'ISO 3310-1 R20/3': ISO_3310_1_R20_3,
                             'ISO 3310-1 R40/3': ISO_3310_1_R40_3,
                             'ISO 3310-1 R10': ISO_3310_1_R10,
                             'ASTM E11': ASTM_E11_sieve_list,}

    def psd_spacing(d_min=None, d_max=None, pts=20, method='logarithmic'):
        r'''Create a particle spacing mesh in one of several ways for use in
        modeling discrete particle size distributions. The allowable meshes are
        'linear', 'logarithmic', a geometric series specified by a Renard number
        such as 'R10', or the meshes available in one of several sieve standards.

        Parameters
        ----------
        d_min : float, optional
            The minimum diameter at which the mesh starts, [m]
        d_max : float, optional
            The maximum diameter at which the mesh ends, [m]
        pts : int, optional
            The number of points to return for the mesh (note this is not respected
            by sieve meshes), [-]
        method : str, optional
            Either 'linear', 'logarithmic', a Renard number like 'R10' or 'R5' or
            'R2.5', or one of the sieve standards 'ISO 3310-1 R40/3',
            'ISO 3310-1 R20', 'ISO 3310-1 R20/3', 'ISO 3310-1', 'ISO 3310-1 R10',
            'ASTM E11', [-]

        Returns
        -------
        ds : list[float]
            The generated mesh diameters, [m]

        Notes
        -----
        Note that when specifying a Renard series, only one of `d_min` or `d_max` can
        be respected! Provide only one of those numbers.

        Note that when specifying a sieve standard the number of points is not
        respected!

        Examples
        --------
        >>> psd_spacing(d_min=5e-5, d_max=5e-4, method='ISO 3310-1 R20/3')
        [6.3e-05, 9e-05, 0.000125, 0.00018, 0.00025, 0.000355, 0.0005]

        References
        ----------
        .. [1] ASTM E11 - 17 - Standard Specification for Woven Wire Test Sieve
           Cloth and Test Sieves.
        .. [2] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
           -- Part 1: Test Sieves of Metal Wire Cloth.
        '''
        if d_min is not None:
            d_min = float(d_min)
        if d_max is not None:
            d_max = float(d_max)
        if method == 'logarithmic':
            return logspace(log10(d_min), log10(d_max), pts)
        elif method == 'linear':
            return linspace(d_min, d_max, pts)
        elif method[0] in ('R', 'r'):
            ratio = 10**(1.0/float(method[1:]))
            if d_min is not None and d_max is not None:
                raise ValueError('For geometric (Renard) series, only '
                                'one of `d_min` and `d_max` should be provided')
            if d_min is not None:
                ds = [d_min]
                for i in range(pts-1):
                    ds.append(ds[-1]*ratio)
                return ds
            elif d_max is not None:
                ds = [d_max]
                for i in range(pts-1):
                    ds.append(ds[-1]/ratio)
                return list(reversed(ds))
        elif method in sieve_spacing_options:
            l = sieve_spacing_options[method]
            ds = []
            for sieve in l:
               if  d_min <= sieve.opening <= d_max:
                   ds.append(sieve.opening)
            return list(reversed(ds))
        else:
            raise ValueError('Method not recognized')

    def pdf_lognormal(d, d_characteristic, s):
        r'''Calculates the probability density function of a lognormal particle
        distribution given a particle diameter `d`, characteristic particle
        diameter `d_characteristic`, and distribution standard deviation `s`.

        .. math::
            q(d) = \frac{1}{ds\sqrt{2\pi}} \exp\left[-0.5\left(\frac{
            \ln(d/d_{characteristic})}{s}\right)^2\right]

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        d_characteristic : float
            Characteristic particle diameter; often D[3, 3] is used for this
            purpose but not always, [m]
        s : float
            Distribution standard deviation, [-]

        Returns
        -------
        pdf : float
            Lognormal probability density function, [-]

        Notes
        -----
        The characteristic diameter can be in terns of number density (denoted
        :math:`q_0(d)`), length density (:math:`q_1(d)`), surface area density
        (:math:`q_2(d)`), or volume density (:math:`q_3(d)`). Volume density is
        most often used. Interconversions among the distributions is possible but
        tricky.

        The standard distribution (i.e. the one used in Scipy) can perform the same
        computation with  `d_characteristic` as the value of `scale`.

        >>> import scipy.stats
        >>> float(scipy.stats.lognorm.pdf(x=1E-4, s=1.1, scale=1E-5))
        405.5420921

        Scipy's calculation is over 300 times slower however, and this expression
        is numerically integrated so speed is required.

        Examples
        --------
        >>> pdf_lognormal(1E-4, 1E-5, 1.1)
        405.5420921

        References
        ----------
        .. [1] ISO 9276-2:2014 - Representation of Results of Particle Size
           Analysis - Part 2: Calculation of Average Particle Sizes/Diameters and
           Moments from Particle Size Distributions.
        '''
        try:
            log_term = log(d/d_characteristic)/s
        except ValueError:
            return 0.0
        return 1./(d*s*ROOT_TWO_PI)*exp(-0.5*log_term*log_term)

    def cdf_lognormal(d, d_characteristic, s):
        r'''Calculates the cumulative distribution function of a lognormal particle
        distribution given a particle diameter `d`, characteristic particle
        diameter `d_characteristic`, and distribution standard deviation `s`.

        .. math::
            Q(d) = 0.5\left(1 + \text{err}\left[\left(\frac{\ln(d/d_c)}{s\sqrt{2}}
            \right)\right]\right)

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        d_characteristic : float
            Characteristic particle diameter; often D[3, 3] is used for this
            purpose but not always, [m]
        s : float
            Distribution standard deviation, [-]

        Returns
        -------
        cdf : float
            Lognormal cumulative density function, [-]

        Notes
        -----
        The characteristic diameter can be in terns of number density (denoted
        :math:`q_0(d)`), length density (:math:`q_1(d)`), surface area density
        (:math:`q_2(d)`), or volume density (:math:`q_3(d)`). Volume density is
        most often used. Interconversions among the distributions is possible but
        tricky.

        The standard distribution (i.e. the one used in Scipy) can perform the same
        computation with  `d_characteristic` as the value of `scale`.

        >>> import scipy.stats
        >>> float(scipy.stats.lognorm.cdf(x=1E-4, s=1.1, scale=1E-5))
        0.9818369875798177

        Scipy's calculation is over 100 times slower however.

        Examples
        --------
        >>> cdf_lognormal(d=1E-4, d_characteristic=1E-5, s=1.1)
        0.9818369875798

        References
        ----------
        .. [1] ISO 9276-2:2014 - Representation of Results of Particle Size
           Analysis - Part 2: Calculation of Average Particle Sizes/Diameters and
           Moments from Particle Size Distributions.
        '''
        try:
            return 0.5*(1.0 + erf((log(d/d_characteristic))/(s*sqrt(2.0))))
        except:
            # math error at cdf = 0 (x going as low as possible)
            return 0.0

    def pdf_lognormal_basis_integral(d, d_characteristic, s, n):
        r'''Calculates the integral of the multiplication of d^n by the lognormal
        pdf, given a particle diameter `d`, characteristic particle
        diameter `d_characteristic`, distribution standard deviation `s`, and
        exponent `n`.

        .. math::
            \int d^n\cdot q(d)\; dd = -\frac{1}{2} \exp\left(\frac{s^2 n^2}{2}
            \right)d^n \left(\frac{d}{d_{characteristic}}\right)^{-n}
            \text{erf}\left[\frac{s^2 n - \ln(d/d_{characteristic})}
            {\sqrt{2} s} \right]

        This is the crucial integral required for interconversion between different
        bases such as number density (denoted :math:`q_0(d)`), length density
        (:math:`q_1(d)`), surface area density (:math:`q_2(d)`), or volume density
        (:math:`q_3(d)`).

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        d_characteristic : float
            Characteristic particle diameter; often D[3, 3] is used for this
            purpose but not always, [m]
        s : float
            Distribution standard deviation, [-]
        n : int
            Exponent of the multiplied n

        Returns
        -------
        pdf_basis_integral : float
            Integral of lognormal pdf multiplied by d^n, [-]

        Notes
        -----
        This integral has been verified numerically. This integral is itself
        integrated, so it is crucial to obtain an analytical form for at least
        this integral.

        Note overflow or zero division issues may occur for very large values of
        `s`, larger than 10. No mathematical limit was able to be obtained with
        a CAS.

        Examples
        --------
        >>> pdf_lognormal_basis_integral(d=1E-4, d_characteristic=1E-5, s=1.1, n=-2)
        56228306549.26362
        '''
        try:
            s2 = s*s
            t0 = exp(s2*n*n*0.5)
            d_ratio = d/d_characteristic
            t1 = (d/(d_ratio))**n
            t2 = erf((s2*n - log(d_ratio))/(sqrt(2.)*s))
            return -0.5*t0*t1*t2
        except (OverflowError, ZeroDivisionError, ValueError):
            return pdf_lognormal_basis_integral(d=1E-80, d_characteristic=d_characteristic, s=s, n=n)

    def pdf_Gates_Gaudin_Schuhman(d, d_characteristic, m):
        r'''Calculates the probability density of a particle
        distribution following the Gates, Gaudin and Schuhman (GGS) model given a
        particle diameter `d`, characteristic (maximum) particle
        diameter `d_characteristic`, and exponent `m`.

        .. math::
            q(d) = \frac{n}{d}\left(\frac{d}{d_{characteristic}}\right)^m
            \text{ if } d < d_{characteristic} \text{ else } 0

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        d_characteristic : float
            Characteristic particle diameter; in this model, it is the largest
            particle size diameter in the distribution, [m]
        m : float
            Particle size distribution exponent, [-]

        Returns
        -------
        pdf : float
            GGS probability density function, [-]

        Notes
        -----
        The characteristic diameter can be in terns of number density (denoted
        :math:`q_0(d)`), length density (:math:`q_1(d)`), surface area density
        (:math:`q_2(d)`), or volume density (:math:`q_3(d)`). Volume density is
        most often used. Interconversions among the distributions is possible but
        tricky.

        Examples
        --------
        >>> pdf_Gates_Gaudin_Schuhman(d=2E-4, d_characteristic=1E-3, m=2.3)
        283.8355768512045

        References
        ----------
        .. [1] Schuhmann, R., 1940. Principles of Comminution, I-Size Distribution
           and Surface Calculations. American Institute of Mining, Metallurgical
           and Petroleum Engineers Technical Publication 1189. Mining Technology,
           volume 4, p. 1-11.
        .. [2] Bayat, Hossein, Mostafa Rastgo, Moharram Mansouri Zadeh, and Harry
           Vereecken. "Particle Size Distribution Models, Their Characteristics and
           Fitting Capability." Journal of Hydrology 529 (October 1, 2015): 872-89.
        '''
        if d <= d_characteristic:
            return m/d*(d/d_characteristic)**m
        else:
            return 0.0

    def cdf_Gates_Gaudin_Schuhman(d, d_characteristic, m):
        r'''Calculates the cumulative distribution function of a particle
        distribution following the Gates, Gaudin and Schuhman (GGS) model given a
        particle diameter `d`, characteristic (maximum) particle
        diameter `d_characteristic`, and exponent `m`.

        .. math::
            Q(d) = \left(\frac{d}{d_{characteristic}}\right)^m \text{ if }
            d < d_{characteristic} \text{ else } 1

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        d_characteristic : float
            Characteristic particle diameter; in this model, it is the largest
            particle size diameter in the distribution, [m]
        m : float
            Particle size distribution exponent, [-]

        Returns
        -------
        cdf : float
            GGS cumulative density function, [-]

        Notes
        -----
        The characteristic diameter can be in terns of number density (denoted
        :math:`q_0(d)`), length density (:math:`q_1(d)`), surface area density
        (:math:`q_2(d)`), or volume density (:math:`q_3(d)`). Volume density is
        most often used. Interconversions among the distributions is possible but
        tricky.

        Examples
        --------
        >>> cdf_Gates_Gaudin_Schuhman(d=2E-4, d_characteristic=1E-3, m=2.3)
        0.024681354508800397

        References
        ----------
        .. [1] Schuhmann, R., 1940. Principles of Comminution, I-Size Distribution
           and Surface Calculations. American Institute of Mining, Metallurgical
           and Petroleum Engineers Technical Publication 1189. Mining Technology,
           volume 4, p. 1-11.
        .. [2] Bayat, Hossein, Mostafa Rastgo, Moharram Mansouri Zadeh, and Harry
           Vereecken. "Particle Size Distribution Models, Their Characteristics and
           Fitting Capability." Journal of Hydrology 529 (October 1, 2015): 872-89.
        '''
        if d <= d_characteristic:
            return (d/d_characteristic)**m
        else:
            return 1.0

    def pdf_Gates_Gaudin_Schuhman_basis_integral(d, d_characteristic, m, n):
        r'''Calculates the integral of the multiplication of d^n by the Gates,
        Gaudin and Schuhman (GGS) model given a particle diameter `d`,
        characteristic (maximum) particle diameter `d_characteristic`, and exponent
        `m`.

        .. math::
            \int d^n\cdot q(d)\; dd =\frac{m}{m+n} d^n \left(\frac{d}
            {d_{characteristic}}\right)^m

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        d_characteristic : float
            Characteristic particle diameter; in this model, it is the largest
            particle size diameter in the distribution, [m]
        m : float
            Particle size distribution exponent, [-]
        n : int
            Exponent of the multiplied n, [-]

        Returns
        -------
        pdf_basis_integral : float
            Integral of Rosin Rammler pdf multiplied by d^n, [-]

        Notes
        -----
        This integral does not have any numerical issues as `d` approaches 0.

        Examples
        --------
        >>> pdf_Gates_Gaudin_Schuhman_basis_integral(d=2E-4, d_characteristic=1E-3, m=2.3, n=-3)
        -10136984887.543015
        '''
        return m/(m+n)*d**n*(d/d_characteristic)**m

    def pdf_Rosin_Rammler(d, k, m):
        r'''Calculates the probability density of a particle
        distribution following the Rosin-Rammler (RR) model given a
        particle diameter `d`, and the two parameters `k` and `m`.

        .. math::
            q(d) = k m d^{(m-1)} \exp(- k d^{m})

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        k : float
            Parameter in the model, [(1/m)^m]
        m : float
            Parameter in the model, [-]

        Returns
        -------
        pdf : float
            RR probability density function, [-]

        Notes
        -----

        Examples
        --------
        >>> pdf_Rosin_Rammler(1E-3, 200, 2)
        0.3999200079994667

        References
        ----------
        .. [1] Rosin, P. "The Laws Governing the Fineness of Powdered Coal." J.
           Inst. Fuel. 7 (1933): 29-36.
        .. [2] Bayat, Hossein, Mostafa Rastgo, Moharram Mansouri Zadeh, and Harry
           Vereecken. "Particle Size Distribution Models, Their Characteristics and
           Fitting Capability." Journal of Hydrology 529 (October 1, 2015): 872-89.
        '''
        return d**(m - 1.0)*k*m*exp(-d**m*k)

    def cdf_Rosin_Rammler(d, k, m):
        r'''Calculates the cumulative distribution function of a particle
        distribution following the Rosin-Rammler (RR) model given a
        particle diameter `d`, and the two parameters `k` and `m`.

        .. math::
            Q(d) = 1 - \exp\left(-k d^m\right)

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        k : float
            Parameter in the model, [(1/m)^m]
        m : float
            Parameter in the model, [-]

        Returns
        -------
        cdf : float
            RR cumulative density function, [-]

        Notes
        -----
        The characteristic diameter can be in terns of number density (denoted
        :math:`q_0(d)`), length density (:math:`q_1(d)`), surface area density
        (:math:`q_2(d)`), or volume density (:math:`q_3(d)`). Volume density is
        most often used. Interconversions among the distributions is possible but
        tricky.

        Examples
        --------
        >>> cdf_Rosin_Rammler(5E-2, 200, 2)
        0.3934693402873667

        References
        ----------
        .. [1] Rosin, P. "The Laws Governing the Fineness of Powdered Coal." J.
           Inst. Fuel. 7 (1933): 29-36.
        .. [2] Bayat, Hossein, Mostafa Rastgo, Moharram Mansouri Zadeh, and Harry
           Vereecken. "Particle Size Distribution Models, Their Characteristics and
           Fitting Capability." Journal of Hydrology 529 (October 1, 2015): 872-89.
        '''
        return 1.0 - exp(-k*d**m)

    def pdf_Rosin_Rammler_basis_integral(d, k, m, n):
        r'''Calculates the integral of the multiplication of d^n by the Rosin
        Rammler (RR) pdf, given a particle diameter `d`, and the two parameters `k`
        and `m`.

        .. math::
            \int d^n\cdot q(d)\; dd =-d^{m+n} k(d^mk)^{-\frac{m+n}{m}}\Gamma
            \left(\frac{m+n}{m}\right)\text{gammaincc}\left[\left(\frac{m+n}{m}
            \right), kd^m\right]

        Parameters
        ----------
        d : float
            Specified particle diameter, [m]
        k : float
            Parameter in the model, [(1/m)^m]
        m : float
            Parameter in the model, [-]
        n : int
            Exponent of the multiplied n, [-]

        Returns
        -------
        pdf_basis_integral : float
            Integral of Rosin Rammler pdf multiplied by d^n, [-]

        Notes
        -----
        This integral was derived using a CAS, and verified numerically.
        The `gammaincc` function is that from scipy.special, and `gamma` from the
        same.

        For very high powers of `n` or `m` when the diameter is very low,
        exceptions may occur.

        Examples
        --------
        >>> "{:g}".format(pdf_Rosin_Rammler_basis_integral(5E-2, 200, 2, 3))
        '-0.000452399'
        '''
        # Also not able to compute the limit for d approaching 0.
        try:
            a = (m + n)/m
            x = d**m*k
            t1 = float(gamma(a))*float(gammaincc(a, x))
            return (-d**(m+n)*k*(d**m*k)**(-a))*t1
        except (OverflowError, ZeroDivisionError) as e:
            if d == 1E-40:
                raise e
            return pdf_Rosin_Rammler_basis_integral(1E-40, k, m, n)

    names = {0: 'Number distribution', 1: 'Length distribution',
             2: 'Area distribution', 3: 'Volume/Mass distribution'}

    def _label_distribution_n(n):  # pragma: no cover
        if n in names:
            return names[n]
        else:
            return f'Order {n!s} distribution'

    _mean_size_docstring = r"""Calculates the mean particle size according to moment-ratio
            notation. This is the more common and often convenient definition.

            .. math::
                \left[\bar D_{p,q} \right]^{(p-q)} = \frac{\sum_i n_i  D_i^p }
                {\sum_i n_i D_i^q}

                \left[\bar D_{p,p} \right] = \exp\left[\frac{\sum_i n_i  D_i^p\ln
                D_i }{\sum_i n_i D_i^p}\right]  \text{, if p = q}

            Note that :math:`n_i` in the above equation is replaceable with
            the fraction of particles in that bin.

            Parameters
            ----------
            p : int
                Power and/or subscript of D moment in the above equations, [-]
            q : int
                Power and/or subscript of D moment in the above equations, [-]

            Returns
            -------
            d_pq : float
                Mean particle size according to the specified p and q, [m]

            Notes
            -----
            The following is a list of common names for specific mean diameters.

            * **D[-3, 0]**: arithmetic harmonic mean volume diameter
            * **D[-2, 1]**: size-weighted harmonic mean volume diameter
            * **D[-1, 2]**: area-weighted harmonic mean volume diameter
            * **D[-2, 0]**: arithmetic harmonic mean area diameter
            * **D[-1, 1]**: size-weighted harmonic mean area diameter
            * **D[-1, 0]**: arithmetic harmonic mean diameter
            * **D[0, 0]**: arithmetic geometric mean diameter
            * **D[1, 1]**: size-weighted geometric mean diameter
            * **D[2, 2]**: area-weighted geometric mean diameter
            * **D[3, 3]**: volume-weighted geometric mean diameter
            * **D[1, 0]**: arithmetic mean diameter
            * **D[2, 1]**: size-weighted mean diameter
            * **D[3, 2]**: area-weighted mean diameter, **Sauter mean diameter**
            * **D[4, 3]**: volume-weighted mean diameter, **De Brouckere diameter**
            * **D[2, 0]**: arithmetic mean area diameter
            * **D[3, 1]**: size-weighted mean area diameter
            * **D[4, 2]**: area-weighted mean area diameter
            * **D[5, 3]**: volume-weighted mean area diameter
            * **D[3, 0]**: arithmetic mean volume diameter
            * **D[4, 1]**: size-weighted mean volume diameter
            * **D[5, 2]**: area-weighted mean volume diameter
            * **D[6, 3]**: volume-weighted mean volume diameter

            This notation was first introduced in [1]_.

            The sum of p and q is called the order of the mean size [3]_.

            .. math::
                \bar D_{p,q}  \equiv \bar D_{q, p}

            Examples
            --------
    %s

            References
            ----------
            .. [1] Mugele, R. A., and H. D. Evans. "Droplet Size Distribution in
               Sprays." Industrial & Engineering Chemistry 43, no. 6 (June 1951):
               1317-24. https://doi.org/10.1021/ie50498a023.
            .. [2] ASTM E799 - 03(2015) - Standard Practice for Determining Data
               Criteria and Processing for Liquid Drop Size Analysis.
            .. [3] ISO 9276-2:2014 - Representation of Results of Particle Size
               Analysis - Part 2: Calculation of Average Particle Sizes/Diameters
               and Moments from Particle Size Distributions.
    """

    _mean_size_iso_docstring =  r"""Calculates the mean particle size according to moment
            notation (ISO). This system is related to the moment-ratio notation
            as follows; see the `mean_size` method for the full formulas.

            .. math::
                \bar x_{p-q, q} \equiv \bar x_{k+r, r}  \equiv \bar D_{p,q}

            Parameters
            ----------
            k : int
                Power and/or subscript of D moment in the above equations, [-]
            r : int
                Power and/or subscript of D moment in the above equations, [-]

            Returns
            -------
            x_kr : float
                Mean particle size according to the specified k and r in the ISO
                series, [m]

            Notes
            -----
            The following is a list of common names for specific mean diameters in
            the ISO naming convention.

            * **x[-3, 0]**: arithmetic harmonic mean volume diameter
            * **x[-3, 1]**: size-weighted harmonic mean volume diameter
            * **x[-3, 2]**: area-weighted harmonic mean volume diameter
            * **x[-2, 0]**: arithmetic harmonic mean area diameter
            * **x[-2, 1]**: size-weighted harmonic mean area diameter
            * **x[-1, 0]**: arithmetic harmonic mean diameter
            * **x[0, 0]**: arithmetic geometric mean diameter
            * **x[0, 1]**: size-weighted geometric mean diameter
            * **x[0, 2]**: area-weighted geometric mean diameter
            * **x[0, 3]**: volume-weighted geometric mean diameter
            * **x[1, 0]**: arithmetic mean diameter
            * **x[1, 1]**: size-weighted mean diameter
            * **x[1, 2]**: area-weighted mean diameter, **Sauter mean diameter**
            * **x[1, 3]**: volume-weighted mean diameter, **De Brouckere diameter**
            * **x[2, 0]**: arithmetic mean area diameter
            * **x[1, 1]**: size-weighted mean area diameter
            * **x[2, 2]**: area-weighted mean area diameter
            * **x[2, 3]**: volume-weighted mean area diameter
            * **x[3, 0]**: arithmetic mean volume diameter
            * **x[3, 1]**: size-weighted mean volume diameter
            * **x[3, 2]**: area-weighted mean volume diameter
            * **x[3, 3]**: volume-weighted mean volume diameter

            When working with continuous distributions, the ISO series must be used
            to perform the actual calculations.

            Examples
            --------
    %s

            References
            ----------
            .. [1] ISO 9276-2:2014 - Representation of Results of Particle Size
               Analysis - Part 2: Calculation of Average Particle Sizes/Diameters
               and Moments from Particle Size Distributions.
            """

    class ParticleSizeDistributionContinuous:
        r'''Base class representing a continuous particle size distribution
        specified by a mathematical/statistical function. This class holds the
        common methods only.

        Notes
        -----
        Although the stated units of input are in meters, this class is actually
        independent of the units provided; all results will be consistent with the
        provided unit.

        Examples
        --------
        Example problem from [1]_.

        >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6)

        References
        ----------
        .. [1] ISO 9276-2:2014 - Representation of Results of Particle Size
           Analysis - Part 2: Calculation of Average Particle Sizes/Diameters and
           Moments from Particle Size Distributions.
        '''

        def _pdf_basis_integral_definite(self, d_min, d_max, n):
            # Needed as an api for numerical integrals
            return (self._pdf_basis_integral(d=d_max, n=n)
                    - self._pdf_basis_integral(d=d_min, n=n))

        def pdf(self, d, n=None):
            r'''Computes the probability density function of a
            continuous particle size distribution at a specified particle diameter,
            an optionally in a specified basis. The evaluation function varies with
            the distribution chosen. The interconversion between distribution
            orders is performed using the following formula [1]_:

            .. math::
                q_s(d) = \frac{x^{(s-r)} q_r(d) dd}
                { \int_0^\infty d^{(s-r)} q_r(d) dd}

            Parameters
            ----------
            d : float
                Particle size diameter, [m]
            n : int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer, [-]

            Returns
            -------
            pdf : float
                The probability density function at the specified diameter and
                order, [-]

            Notes
            -----
            The pdf order conversions are typically available analytically after
            some work. They have been verified numerically. See the various
            functions with names ending with 'basis_integral' for the formulations.
            The distributions normally do not have analytical limits for diameters
            of 0 or infinity, but large values suffice to capture the area of the
            integral.

            Examples
            --------
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6, order=3)
            >>> psd.pdf(1e-5)
            30522.765209509154
            >>> psd.pdf(1e-5, n=3)
            30522.765209509154
            >>> psd.pdf(1e-5, n=0)
            1238.661379483343

            References
            ----------
            .. [1] Masuda, Hiroaki, Ko Higashitani, and Hideto Yoshida. Powder
               Technology: Fundamentals of Particles, Powder Beds, and Particle
               Generation. CRC Press, 2006.
            '''
            ans = self._pdf(d=d)
            if n is not None and n != self.order:
                power = n - self.order
                numerator = d**power*ans
                denominator = self._pdf_basis_integral_definite(d_min=0.0, d_max=self.d_excessive, n=power)
                ans = numerator/denominator
            # Handle splines which might go below zero
            ans = max(ans, 0.0)
            if self.truncated:
                if d < self.d_min or d > self.d_max:
                    return 0.0
                ans = (ans)/(self._cdf_d_max - self._cdf_d_min)
            return ans

        def cdf(self, d, n=None):
            r'''Computes the cumulative distribution density function of a
            continuous particle size distribution at a specified particle diameter,
            an optionally in a specified basis. The evaluation function varies with
            the distribution chosen.

            .. math::
                Q_n(d) = \int_0^d q_n(d)

            Parameters
            ----------
            d : float
                Particle size diameter, [m]
            n : int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer, [-]

            Returns
            -------
            cdf : float
                The cumulative distribution function at the specified diameter and
                order, [-]

            Notes
            -----
            Analytical integrals can be found for most distributions even when
            order conversions are necessary.

            Examples
            --------
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6, order=3)
            >>> [psd.cdf(5e-6, n) for n in range(4)]
            [0.933192798731, 0.8413447460685, 0.6914624612740, 0.5]
            '''
            if n is not None and n != self.order:
                power = n - self.order
                # One of the pdf_basis_integral calls could be saved except for
                # support for numerical integrals
                numerator = self._pdf_basis_integral_definite(d_min=0.0, d_max=d, n=power)
                denominator = self._pdf_basis_integral_definite(d_min=0.0, d_max=self.d_excessive, n=power)
                ans =  max(numerator/denominator, 0.0)
            # Handle splines which might go below zero
            else:
                ans = max(self._cdf(d=d), 0.0)
            if self.truncated:
                if d <= self.d_min:
                    return 0.0
                elif d >= self.d_max:
                    return 1.0
                ans = (ans - self._cdf_d_min)/(self._cdf_d_max - self._cdf_d_min)
            return ans

        def delta_cdf(self, d_min, d_max, n=None):
            r'''Computes the difference in cumulative distribution function between
            two particle size diameters.

            .. math::
                \Delta Q_n = Q_n(d_{max}) - Q_n(d_{min})

            Parameters
            ----------
            d_min : float
                Lower particle size diameter, [m]
            d_max : float
                Upper particle size diameter, [m]
            n : int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer, [-]

            Returns
            -------
            delta_cdf : float
                The difference in the cumulative distribution function for the two
                diameters specified, [-]

            Examples
            --------
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6, order=3)
            >>> psd.delta_cdf(1e-6, 1e-5)
            0.9165280099853876
            '''
            return self.cdf(d_max, n=n) - self.cdf(d_min, n=n)

        def dn(self, fraction, n=None):
            r'''Computes the diameter at which a specified `fraction` of the
            distribution falls under. Utilizes a bounded solver to search for the
            desired diameter.

            Parameters
            ----------
            fraction : float
                Fraction of the distribution which should be under the calculated
                diameter, [-]
            n : int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer, [-]

            Returns
            -------
            d : float
                Particle size diameter, [m]

            Examples
            --------
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6, order=3)
            >>> psd.dn(.5)
            5e-06
            >>> psd.dn(1)
            0.0002947436533523378
            >>> psd.dn(0)
            0.0
            '''
            fraction = float(fraction)
            if fraction == 1.0:
                # Avoid returning the maximum value of the search interval
                fraction = 1.0 - epsilon
            if fraction < 0:
                raise ValueError('Fraction must be more than 0')
            elif fraction == 0:  # pragma: no cover
                if self.truncated:
                    return self.d_min
                return 0.0
                # Solve to float prevision limit - works well, but is there a real
                # point when with mpmath it would never happen?
                # dist.cdf(dist.dn(0)-1e-35) == 0
                # dist.cdf(dist.dn(0)-1e-36) == input
                # dn(0) == 1.9663615597466143e-20
    #            def err(d):
    #                cdf = self.cdf(d, n=n)
    #                if cdf == 0:
    #                    cdf = -1
    #                return cdf
    #            return brenth(err, self.d_minimum, self.d_excessive, maxiter=1000, xtol=1E-200)

            elif fraction > 1:
                raise ValueError('Fraction less than 1')
            # As the dn may be incredibly small, it is required for the absolute
            # tolerance to not be happy - it needs to continue iterating as long
            # as necessary to pin down the answer
            return brenth(lambda d:self.cdf(d, n=n) -fraction,
                          self.d_minimum, self.d_excessive, maxiter=1000, xtol=1E-200)

        def ds_discrete(self, d_min=None, d_max=None, pts=20, limit=1e-9,
                        method='logarithmic'):
            r'''Create a particle spacing mesh to perform calculations with,
            according to one of several ways. The allowable meshes are
            'linear', 'logarithmic', a geometric series specified by a Renard
            number such as 'R10', or the meshes available in one of several sieve
            standards.

            Parameters
            ----------
            d_min : float, optional
                The minimum diameter at which the mesh starts, [m]
            d_max : float, optional
                The maximum diameter at which the mesh ends, [m]
            pts : int, optional
                The number of points to return for the mesh (note this is not
                respected by sieve meshes), [-]
            limit : float
                If `d_min` or `d_max` is not specified, it will be calculated as the
                `dn` at which this limit or 1-limit exists (this is ignored for
                Renard numbers), [-]
            method : str, optional
                Either 'linear', 'logarithmic', a Renard number like 'R10' or 'R5'
                or'R2.5', or one of the sieve standards 'ISO 3310-1 R40/3',
                'ISO 3310-1 R20', 'ISO 3310-1 R20/3', 'ISO 3310-1',
                'ISO 3310-1 R10', 'ASTM E11', [-]

            Returns
            -------
            ds : list[float]
                The generated mesh diameters, [m]

            Notes
            -----
            Note that when specifying a Renard series, only one of `d_min` or `d_max` can
            be respected! Provide only one of those numbers.

            Note that when specifying a sieve standard the number of points is not
            respected!

            References
            ----------
            .. [1] ASTM E11 - 17 - Standard Specification for Woven Wire Test Sieve
               Cloth and Test Sieves.
            .. [2] ISO 3310-1:2016 - Test Sieves -- Technical Requirements and Testing
               -- Part 1: Test Sieves of Metal Wire Cloth.
            '''
            if method[0] not in ('R', 'r'):
                if d_min is None:
                    d_min = self.dn(limit)
                if d_max is None:
                    d_max = self.dn(1.0 - limit)
            return psd_spacing(d_min=d_min, d_max=d_max, pts=pts, method=method)

        def fractions_discrete(self, ds, n=None):
            r'''Computes the fractions of the cumulative distribution functions
            which lie between the specified specified particle diameters. The first
            diameter contains the cdf from 0 to it.

            Parameters
            ----------
            ds : list[float]
                Particle size diameters, [m]
            n : int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer, [-]

            Returns
            -------
            fractions : float
                The differences in the cumulative distribution functions at the
                specified diameters and order, [-]

            Examples
            --------
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6, order=3)
            >>> psd.fractions_discrete([1e-6, 1e-5, 1e-4, 1e-3])
            [0.00064347101291, 0.916528009985, 0.0828285179619, 1.039798e-09]
            '''
            cdfs = [self.cdf(d, n=n) for d in ds]
            return [cdfs[0]] + diff(cdfs)

        def cdf_discrete(self, ds, n=None):
            r'''Computes the cumulative distribution functions for a list of
            specified particle diameters.

            Parameters
            ----------
            ds : list[float]
                Particle size diameters, [m]
            n : int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer, [-]

            Returns
            -------
            cdfs : float
                The cumulative distribution functions at the specified diameters
                and order, [-]

            Examples
            --------
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6, order=3)
            >>> psd.cdf_discrete([1e-6, 1e-5, 1e-4, 1e-3])
            [0.000643471012913, 0.917171480998, 0.999999998960, 1.0]
            '''
            return [self.cdf(d, n=n) for d in ds]

        def mean_size(self, p, q):
            '''
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6)
            >>> psd.mean_size(3, 2)
            4.412484512922977e-06

            Note that for the case where p == q, a different set of formulas are
            required - which do not have analytical results for many distributions.
            Therefore, a close numerical approximation is used instead, to
            perturb the values of p and q so they are 1E-9 away from each other.
            This leads only to slight errors, as in the example below where the
            correct answer is 5E-6.

            >>> psd.mean_size(3, 3)
            4.9999999304923345e-06
            '''
            if p == q:
                p -= 1e-9
                q += 1e-9
            pow1 = q - self.order

            denominator = self._pdf_basis_integral_definite(d_min=self.d_minimum, d_max=self.d_excessive, n=pow1)
            root_power = p - q
            pow3 = p - self.order
            numerator = self._pdf_basis_integral_definite(d_min=self.d_minimum, d_max=self.d_excessive, n=pow3)
            return float((numerator/denominator)**(1.0/(root_power)))

        def mean_size_ISO(self, k, r):
            '''
            >>> psd = PSDLognormal(s=0.5, d_characteristic=5E-6)
            >>> psd.mean_size_ISO(1, 2)
            4.412484512922977e-06
            '''
            p = k + r
            q = r
            return self.mean_size(p=p, q=q)

        @property
        def vssa(self):
            r'''The volume-specific surface area of a particle size distribution.

            .. math::
                \text{VSSA} = \frac{6}{\bar x_{1,2}}

            Returns
            -------
            VSSA : float
                The volume-specific surface area of the distribution, [m^2/m^3]

            Examples
            --------
            >>> PSDLognormal(s=0.5, d_characteristic=5E-6).vssa
            1359778.1436801916

            References
            ----------
            .. [1] ISO 9276-2:2014 - Representation of Results of Particle Size
               Analysis - Part 2: Calculation of Average Particle Sizes/Diameters
               and Moments from Particle Size Distributions.
            '''
            return 6/self.mean_size(3, 2)

        def plot_pdf(self, n=(0, 1, 2, 3), d_min=None, d_max=None, pts=500,
                     normalized=False, method='linear'):  # pragma: no cover
            r'''Plot the probability density function of the particle size
            distribution. The plotted range can be specified using `d_min` and
            `d_max`, or estimated automatically. One or more order can be plotted,
            by providing an iterable of ints as the value of `n` or just one int.

            Parameters
            ----------
            n : tuple(int) or int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer; as many as desired may be specified, [-]
            d_min : float, optional
                Lower particle size diameter, [m]
            d_max : float, optional
                Upper particle size diameter, [m]
            pts : int, optional
                The number of points for values to be calculated, [-]
            normalized : bool, optional
                Whether to display the actual probability density function, which
                may have a huge magnitude - or to divide each point by the sum
                of all the points. Doing this is a common practice, but the values
                at each point are dependent on the number of points being plotted,
                and the distribution of the points;
                [-]
            method : str, optional
                Either 'linear', 'logarithmic', a Renard number like 'R10' or 'R5'
                or'R2.5', or one of the sieve standards 'ISO 3310-1 R40/3',
                'ISO 3310-1 R20', 'ISO 3310-1 R20/3', 'ISO 3310-1',
                'ISO 3310-1 R10', 'ASTM E11', [-]
            '''
            try:
                import matplotlib.pyplot as plt
            except:  # pragma: no cover
                raise ValueError(NO_MATPLOTLIB_MSG)
            ds = self.ds_discrete(d_min=d_min, d_max=d_max, pts=pts, method=method)
            try:
                for ni in n:
                    fractions = [self.pdf(d, n=ni) for d in ds]
                    if normalized:
                        fractions = normalize(fractions)
                    plt.semilogx(ds, fractions, label=_label_distribution_n(ni))
            except Exception:
                fractions = [self.pdf(d, n=n) for d in ds]
                if normalized:
                    fractions = normalize(fractions)
                plt.semilogx(ds, fractions, label=_label_distribution_n(n))
            plt.ylabel('Probability density function, [-]')
            plt.xlabel('Particle diameter, [m]')
            plt.title(f'Probability density function of {self.name} distribution with '
                      f'parameters {self.parameters}')
            plt.legend()
            plt.show()
            return fractions

        def plot_cdf(self, n=(0, 1, 2, 3), d_min=None, d_max=None, pts=500,
                     method='logarithmic'):   # pragma: no cover
            r'''Plot the cumulative distribution function of the particle size
            distribution. The plotted range can be specified using `d_min` and
            `d_max`, or estimated automatically. One or more order can be plotted,
            by providing an iterable of ints as the value of `n` or just one int.

            Parameters
            ----------
            n : tuple(int) or int, optional
                None (for the `order` specified when the distribution was created),
                0 (number), 1 (length), 2 (area), 3 (volume/mass),
                or any integer; as many as desired may be specified, [-]
            d_min : float, optional
                Lower particle size diameter, [m]
            d_max : float, optional
                Upper particle size diameter, [m]
            pts : int, optional
                The number of points for values to be calculated, [-]
            method : str, optional
                Either 'linear', 'logarithmic', a Renard number like 'R10' or 'R5'
                or'R2.5', or one of the sieve standards 'ISO 3310-1 R40/3',
                'ISO 3310-1 R20', 'ISO 3310-1 R20/3', 'ISO 3310-1',
                'ISO 3310-1 R10', 'ASTM E11', [-]
            '''
            try:
                import matplotlib.pyplot as plt
            except:  # pragma: no cover
                raise ValueError(NO_MATPLOTLIB_MSG)

            ds = self.ds_discrete(d_min=d_min, d_max=d_max, pts=pts, method=method)
            try:
                for ni in n:
                    cdfs = self.cdf_discrete(ds=ds, n=ni)
                    plt.semilogx(ds, cdfs, label=_label_distribution_n(ni))
            except:
                cdfs = self.cdf_discrete(ds=ds, n=n)
                plt.semilogx(ds, cdfs, label=_label_distribution_n(n))
            if self.points:
                plt.plot(self.ds, self.fraction_cdf, '+', label='Volume/Mass points')

                if hasattr(self, 'area_fractions'):
                    plt.plot(self.ds, cumsum(self.area_fractions), '+', label='Area points')
                if hasattr(self, 'length_fractions'):
                    plt.plot(self.ds, cumsum(self.length_fractions), '+', label='Length points')
                if hasattr(self, 'number_fractions'):
                    plt.plot(self.ds, cumsum(self.number_fractions), '+', label='Number points')

            plt.ylabel('Cumulative density function, [-]')
            plt.xlabel('Particle diameter, [m]')
            plt.title(f'Cumulative density function of {self.name} distribution with '
                      f'parameters {self.parameters}')
            plt.legend()
            plt.show()

    class ParticleSizeDistribution(ParticleSizeDistributionContinuous):
        r'''Class representing a discrete particle size distribution specified by a
        series of diameter bins, and the quantity of particles in each bin. The
        quantities may be specified as either the fraction of particles in each
        bin, or as cumulative distributions. The input fractions can be
        specified to be in a mass basis (`order=3`), number basis (`order=0`),
        or the orders in between for length basis or area basis. If the
        fractions do not sum to 1, and `cdf` is False, then the fractions are
        normalized. This allows flow rates or counts of size bins to be given as
        well.

        Parameters
        ----------
        ds : list[float]
            Diameter bins; length of the specified quantities, optionally +1 that
            length to specify a cutoff diameter for the smallest diameter bin, [m]
        fractions : list[float], optional
            The mass/mole/volume/length/area/count fractions or cumulative
            distributions or counts of each particle size in
            each diameter bin (the type is specified by `order`), [-]
        order : int, optional
            0 for a number distribution as input; 1 for length distribution;
            2 for area distribution; 3 for mass, mole, or volume distribution, [-]
        cdf : bool, optional
            If the distribution is given as increasing fractions with 1 as the last
            result, `cdf` must be set to True, [-]
        monotonic : bool, optional
            If True, for interpolated quanties, monotonic splines will be used
            instead of the standard splines, [-]

        Attributes
        ----------
        fractions : list[float]
            The mass/mole/volume basis fractions of particles in each bin, [-]
        area_fractions : list[float]
            The area fractions of particles in each bin, [-]
        length_fractions : list[float]
            The length fractions of particles in each bin, [-]
        number_fractions : list[float]
            The number fractions of particles in each bin, [-]
        fraction_cdf : list[float]
            The cumulative mass/mole/volume basis fractions of particles in each
            bin, [-]
        area_cdf : list[float]
            The cumulative area fractions of particles in each bin, [-]
        length_cdf : list[float]
            The cumulative length fractions of particles in each bin, [-]
        number_cdf : list[float]
            The cumulative number fractions of particles in each bin, [-]
        size_classes : bool
            Whether or not the diameter bins were set as size classes (as length
            of fractions + 1), [-]
        N : int
            The number of provided points, [-]

        Notes
        -----
        Although the stated units of input are in meters, this class is actually
        independent of the units provided; all results will be consistent with the
        provided unit.

        Examples
        --------
        Example problem from [1]_, calculating several diameters and the cumulative
        distribution.

        >>> import numpy as np
        >>> ds = 1E-6*np.array([240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532])
        >>> numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]
        >>> psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0)
        >>> psd
        <Particle Size Distribution, points=14, D[3, 3]=0.002451 m>

        References
        ----------
        .. [1] ASTM E799 - 03(2015) - Standard Practice for Determining Data
           Criteria and Processing for Liquid Drop Size Analysis.
        .. [2] ISO 9276-2:2014 - Representation of Results of Particle Size
           Analysis - Part 2: Calculation of Average Particle Sizes/Diameters and
           Moments from Particle Size Distributions.
        '''

        def __repr__(self):
            txt = '<Particle Size Distribution, points=%d, D[3, 3]=%f m>'
            return txt %(self.N, self.mean_size(p=3, q=3))

        size_classes = False
        _interpolated = None
        points = True
        truncated = False
        name = 'Discrete'
        def __init__(self, ds, fractions, cdf=False, order=3, monotonic=True):
            self.monotonic = monotonic
            self.ds = ds
            self.order = order

            if ds is not None and (len(ds) == len(fractions) + 1):
                self.size_classes = True
            else:
                self.size_classes = False

            if cdf:
                # Convert a cdf to fraction set
                if len(fractions)+1 == len(ds):
                    fractions = [fractions[0]] + diff(fractions)
                else:
                    fractions = diff(fractions)
                    fractions.insert(0, 0.0)
            elif sum(fractions) != 1.0:
                # Normalize flow inputs
                tot_inv = 1.0/sum(fractions)
                fractions = [i*tot_inv if i != 0.0 else 0.0 for i in fractions]

            self.N = len(fractions)

            # This will always be in base-3 basis
            if self.order != 3:
                power = 3 - self.order
                d3s = [self.di_power(i, power=power)*fractions[i] for i in range(self.N)]
                tot_d3 = sum(d3s)
                self.fractions = [i/tot_d3 for i in d3s]
            else:
                self.fractions = fractions
            # Set the number fractions
            D3s = [self.di_power(i, power=3) for i in range(self.N)]
            numbers = [Vi/Vp for Vi, Vp in zip(self.fractions, D3s)]
            number_sum = sum(numbers)
            self.number_fractions = [i/number_sum for i in numbers]

            # Set the length fractions
            D3s = [self.di_power(i, power=2) for i in range(self.N)]
            numbers = [Vi/Vp for Vi, Vp in zip(self.fractions, D3s)]
            number_sum = sum(numbers)
            self.length_fractions = [i/number_sum for i in numbers]

            # Set the surface area fractions
            D3s = [self.di_power(i, power=1) for i in range(self.N)]
            numbers = [Vi/Vp for Vi, Vp in zip(self.fractions, D3s)]
            number_sum = sum(numbers)
            self.area_fractions = [i/number_sum for i in numbers]

            # Things for interoperability with the Continuous distribution
            self.d_excessive = float(self.ds[-1])
            self.d_minimum = 0.0
            self.parameters = {}
            self.order = 3
            self.fraction_cdf = self.volume_cdf = cumsum(self.fractions)
            self.area_cdf = cumsum(self.area_fractions)
            self.length_cdf = cumsum(self.length_fractions)
            self.number_cdf = cumsum(self.number_fractions)

        @property
        def interpolated(self):
            if not self._interpolated:
                self._interpolated = PSDInterpolated(ds=self.ds,
                                                     fractions=self.fractions,
                                                     order=3,
                                                     monotonic=self.monotonic)
            return self._interpolated

        def _pdf(self, d):
            return self.interpolated._pdf(d)

        def _cdf(self, d):
            return self.interpolated._cdf(d)

        def _pdf_basis_integral(self, d, n):
            return self.interpolated._pdf_basis_integral(d, n)

        def _fit_obj_function(self, vals, distribution, n):
            err = 0.0
            dist = distribution(*list(vals))
            l = len(self.fractions) if self.size_classes else len(self.fractions) - 1
            for i in range(l):
                delta_cdf = dist.delta_cdf(d_min=self.ds[i], d_max=self.ds[i+1])
                err += abs(delta_cdf - self.fractions[i])
            return err

        def fit(self, x0=None, distribution='lognormal', n=None, **kwargs):
            """Incomplete method to fit experimental values to a curve.

            It is very hard to get good initial guesses, which are really required
            for this. Differential evolution is promising. This API is likely to
            change in the future.
            """
            dist = {'lognormal': PSDLognormal,
                    'GGS': PSDGatesGaudinSchuhman,
                    'RR': PSDRosinRammler}[distribution]

            if distribution == 'lognormal':
                if x0 is None:
                    d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                    s = 0.4
                    x0 = [d_characteristic, s]
            elif distribution == 'GGS':
                if x0 is None:
                    d_characteristic = sum([fi*di for fi, di in zip(self.fractions, self.Dis)])
                    m = 1.5
                    x0 = [d_characteristic, m]
            elif distribution == 'RR' and x0 is None:
                x0 = [5E-6, 1e-2]
            from scipy.optimize import minimize
            return minimize(self._fit_obj_function, x0, args=(dist, n), **kwargs)

        @property
        def Dis(self):
            """Representative diameters of each bin."""
            return [self.di_power(i, power=1) for i in range(self.N)]

        def di_power(self, i, power=1):
            r'''Method to calculate a power of a particle class/bin in a generic
            way so as to support when there are as many `ds` as `fractions`,
            or one more diameter spec than `fractions`.

            When each bin has a lower and upper bound, the formula is as follows
            [1]_.

            .. math::
                D_i^r = \frac{D_{i, ub}^{(r+1)} - D_{i, lb}^{(r+1)}}
                {(D_{i, ub} - D_{i, lb})(r+1)}

            Where `ub` represents the upper bound, and `lb` represents the lower
            bound. Otherwise, the standard definition is used:

            .. math::
                D_i^r = D_i^r

            Parameters
            ----------
            i : int
                The index of the diameter for the calculation, [-]
            power : int
                The exponent, [-]

            Returns
            -------
            di_power : float
                The representative bin diameter raised to  `power`, [m^power]

            References
            ----------
            .. [1] ASTM E799 - 03(2015) - Standard Practice for Determining Data
               Criteria and Processing for Liquid Drop Size Analysis.
            '''
            if self.size_classes:
                rt = power + 1
                return ((self.ds[i+1]**rt - self.ds[i]**rt)/((self.ds[i+1] - self.ds[i])*rt))
            else:
                return self.ds[i]**power

        def mean_size(self, p, q):
            '''
            >>> import numpy as np
            >>> ds = 1E-6*np.array([240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532])
            >>> numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]
            >>> psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0)
            >>> psd.mean_size(3, 2)
            0.002269321031745045
            '''
            if p != q:
                # Note: D(p, q) = D(q, p); in ISO and proven experimentally
                numerator = sum(self.di_power(i=i, power=p)*self.number_fractions[i] for i in range(self.N))
                denominator = sum(self.di_power(i=i, power=q)*self.number_fractions[i] for i in range(self.N))
                return float((numerator/denominator)**(1.0/(p-q)))
            else:
                numerator = sum(log(self.di_power(i=i, power=1))*self.di_power(i=i, power=p)*self.number_fractions[i] for i in range(self.N))
                denominator = sum(self.di_power(i=i, power=q)*self.number_fractions[i] for i in range(self.N))
                return exp(numerator/denominator)

        def mean_size_ISO(self, k, r):
            r'''
            >>> import numpy as np
            >>> ds = 1E-6*np.array([240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532])
            >>> numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]
            >>> psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0)
            >>> psd.mean_size_ISO(1, 2)
            0.002269321031745045
            '''
            p = k + r
            q = r
            return self.mean_size(p=p, q=q)

        @property
        def vssa(self):
            r'''The volume-specific surface area of a particle size distribution.
            Note this uses the diameters provided by the method `Dis`.

            .. math::
                \text{VSSA} = \sum_i \text{fraction}_i \frac{SA_i}{V_i}

            Returns
            -------
            VSSA : float
                The volume-specific surface area of the distribution, [m^2/m^3]

            References
            ----------
            .. [1] ISO 9276-2:2014 - Representation of Results of Particle Size
               Analysis - Part 2: Calculation of Average Particle Sizes/Diameters
               and Moments from Particle Size Distributions.
            '''
            ds = self.Dis
            Vs = [pi/6*di**3 for di in ds]
            SAs = [pi*di**2 for di in ds]
            SASs = [SA/V for SA, V in zip(SAs, Vs)]
            VSSA = sum([fi*SASi for fi, SASi in zip(self.fractions, SASs)])
            return VSSA

    try:  # pragma: no cover
        # Python 2
        ParticleSizeDistributionContinuous.mean_size.__func__.__doc__ = _mean_size_docstring %(ParticleSizeDistributionContinuous.mean_size.__func__.__doc__)
        ParticleSizeDistributionContinuous.mean_size_ISO.__func__.__doc__ = _mean_size_iso_docstring %(ParticleSizeDistributionContinuous.mean_size_ISO.__func__.__doc__)
        ParticleSizeDistribution.mean_size.__func__.__doc__ = _mean_size_docstring %(ParticleSizeDistribution.mean_size.__func__.__doc__)
        ParticleSizeDistribution.mean_size_ISO.__func__.__doc__ = _mean_size_iso_docstring %(ParticleSizeDistribution.mean_size_ISO.__func__.__doc__)
    except AttributeError:  # pragma: no cover
        try:
            # Python 3
            ParticleSizeDistributionContinuous.mean_size.__doc__ = _mean_size_docstring %(ParticleSizeDistributionContinuous.mean_size.__doc__)
            ParticleSizeDistributionContinuous.mean_size_ISO.__doc__ = _mean_size_iso_docstring %(ParticleSizeDistributionContinuous.mean_size_ISO.__doc__)
            ParticleSizeDistribution.mean_size.__doc__ = _mean_size_docstring %(ParticleSizeDistribution.mean_size.__doc__)
            ParticleSizeDistribution.mean_size_ISO.__doc__ = _mean_size_iso_docstring %(ParticleSizeDistribution.mean_size_ISO.__doc__)
        except:
            pass # micropython
    del _mean_size_iso_docstring
    del _mean_size_docstring

    class PSDLognormal(ParticleSizeDistributionContinuous):
        name = 'Lognormal'
        points = False
        truncated = False
        def __init__(self, d_characteristic, s, order=3, d_min=None, d_max=None):
            self.s = s
            self.d_characteristic = d_characteristic
            self.order = order
            self.parameters = {'s': s, 'd_characteristic': d_characteristic,
                               'd_min': d_min, 'd_max': d_max}
            self.d_min = d_min
            self.d_max = d_max
            # Pick an upper bound for the search algorithm of 15 orders of magnitude larger than
            # the characteristic diameter; should never be a problem, as diameters can only range
            # so much, physically.
            if self.d_max is not None:
                self.d_excessive = self.d_max
            else:
                self.d_excessive = 1E15*self.d_characteristic
            if self.d_min is not None:
                self.d_minimum = self.d_min
            else:
                self.d_minimum = 0.0

            if self.d_min is not None or self.d_max is not None:
                self.truncated = True
                if self.d_max is None:
                    self.d_max = self.d_excessive
                if self.d_min is None:
                    self.d_min = 0.0

                self._cdf_d_max = self._cdf(self.d_max)
                self._cdf_d_min = self._cdf(self.d_min)

        def _pdf(self, d):
            return pdf_lognormal(d, self.d_characteristic, self.s)

        def _cdf(self, d):
            return cdf_lognormal(d, self.d_characteristic, self.s)

        def _pdf_basis_integral(self, d, n):
            return pdf_lognormal_basis_integral(d, self.d_characteristic, self.s, n)

    class PSDGatesGaudinSchuhman(ParticleSizeDistributionContinuous):
        name = 'Gates Gaudin Schuhman'
        points = False
        truncated = False
        def __init__(self, d_characteristic, m, order=3, d_min=None, d_max=None):
            self.m = m
            self.d_characteristic = d_characteristic
            self.order = order
            self.parameters = {'m': m, 'd_characteristic': d_characteristic,
                               'd_min': d_min, 'd_max': d_max}

            if self.d_max is not None:
                # PDF above this is zero
                self.d_excessive = self.d_max
            else:
                self.d_excessive = self.d_characteristic
            if self.d_min is not None:
                self.d_minimum = self.d_min
            else:
                self.d_minimum = 0.0

            if self.d_min is not None or self.d_max is not None:
                self.truncated = True
                if self.d_max is None:
                    self.d_max = self.d_excessive
                if self.d_min is None:
                    self.d_min = 0.0

                self._cdf_d_max = self._cdf(self.d_max)
                self._cdf_d_min = self._cdf(self.d_min)

        def _pdf(self, d):
            return pdf_Gates_Gaudin_Schuhman(d, d_characteristic=self.d_characteristic, m=self.m)

        def _cdf(self, d):
            return cdf_Gates_Gaudin_Schuhman(d, d_characteristic=self.d_characteristic, m=self.m)

        def _pdf_basis_integral(self, d, n):
            return pdf_Gates_Gaudin_Schuhman_basis_integral(d, d_characteristic=self.d_characteristic, m=self.m, n=n)

    class PSDRosinRammler(ParticleSizeDistributionContinuous):
        name = 'Rosin Rammler'
        points = False
        truncated = False
        def __init__(self, k, m, order=3, d_min=None, d_max=None):
            self.m = m
            self.k = k
            self.order = order
            self.parameters = {'m': m, 'k': k, 'd_min': d_min, 'd_max': d_max}

            if self.d_max is not None:
                self.d_excessive = self.d_max
            else:
                self.d_excessive = 1e15
            if self.d_min is not None:
                self.d_minimum = self.d_min
            else:
                self.d_minimum = 0.0

            if self.d_min is not None or self.d_max is not None:
                self.truncated = True
                if self.d_max is None:
                    self.d_max = self.d_excessive
                if self.d_min is None:
                    self.d_min = 0.0

                self._cdf_d_max = self._cdf(self.d_max)
                self._cdf_d_min = self._cdf(self.d_min)

        def _pdf(self, d):
            return pdf_Rosin_Rammler(d, k=self.k, m=self.m)

        def _cdf(self, d):
            return cdf_Rosin_Rammler(d, k=self.k, m=self.m)

        def _pdf_basis_integral(self, d, n):
            return pdf_Rosin_Rammler_basis_integral(d, k=self.k, m=self.m, n=n)

    """# These are all brutally slow!
    from scipy.stats import *
    from fluids import *
    distribution = lognorm(s=0.5, scale=5E-6)
    psd = PSDCustom(distribution)

    # psd.dn(0.5, n=2.0) # Doesn't work at all, but the main things do including plots
    """

    class PSDCustom(ParticleSizeDistributionContinuous):
        name = ''
        points = False
        truncated = False
        def __init__(self, distribution, order=3.0, d_excessive=1.0, name=None,
                     d_min=None, d_max=None):
            if name:
                self.name = name
            else:
                try:
                    self.name = distribution.dist.__class__.__name__
                except:
                    pass
            try:
                self.parameters = dict(distribution.kwds)
                self.parameters.update({'d_min': d_min, 'd_max': d_max})
            except:
                self.parameters = {'d_min': d_min, 'd_max': d_max}

            self.distribution = distribution
            self.order = order
            self.d_max = d_max
            self.d_min = d_min

            if self.d_max is not None:
                self.d_excessive = self.d_max
            else:
                self.d_excessive = d_excessive
            if self.d_min is not None:
                self.d_minimum = self.d_min
            else:
                self.d_minimum = 0.0

            if self.d_min is not None or self.d_max is not None:
                self.truncated = True
                if self.d_max is None:
                    self.d_max = self.d_excessive
                if self.d_min is None:
                    self.d_min = 0.0

                self._cdf_d_max = self._cdf(self.d_max)
                self._cdf_d_min = self._cdf(self.d_min)

        def _pdf(self, d):
            return self.distribution.pdf(d)

        def _cdf(self, d):
            return self.distribution.cdf(d)

        def _pdf_basis_integral_definite(self, d_min, d_max, n):
            # Needed as an api for numerical integrals
            n = float(n)
            if d_min == 0:
                d_min = d_max*1E-12

            if n == 0:
                to_int = lambda d : self._pdf(d)
            elif n == 1:
                to_int = lambda d : d*self._pdf(d)
            elif n == 2:
                to_int = lambda d : d*d*self._pdf(d)
            elif n == 3:
                to_int = lambda d : d*d*d*self._pdf(d)
            else:
                to_int = lambda d : d**n*self._pdf(d)

            # points = logspace(log10(max(d_max*1e-3, d_min)), log10(d_max*.999), 40)
            points = [d_max*1e-3] # d_min*.999 d_min
            return float(quad(to_int, d_min, d_max, points=points, epsrel=1e-11)[0])

    class PSDInterpolated(ParticleSizeDistributionContinuous):
        name = 'Interpolated'
        points = True
        truncated = False
        def __init__(self, ds, fractions, order=3, monotonic=True):
            self.order = order
            self.monotonic = monotonic
            self.parameters = {}

            ds = list(ds)
            fractions = list(fractions)

            if len(ds) == len(fractions)+1:
                # size classes, the last point will be zero
                fractions.insert(0, 0.0)
                self.d_minimum = min(ds)
            elif ds[0] != 0:
                ds = [0] + ds
                if len(ds) != len(fractions):
                    fractions = [0] + fractions
                self.d_minimum = 0.0

            self.ds = ds
            self.fractions = fractions

            self.d_excessive = max(ds)

            self.fraction_cdf = cumsum(fractions)
            if self.monotonic:
                from scipy.interpolate import PchipInterpolator
                globals()['PchipInterpolator'] = PchipInterpolator

                self.cdf_spline = PchipInterpolator(ds, self.fraction_cdf, extrapolate=True)
                self.pdf_spline = PchipInterpolator(ds, self.fraction_cdf, extrapolate=True).derivative(1)
            else:
                from scipy.interpolate import UnivariateSpline
                globals()['UnivariateSpline'] = UnivariateSpline

                self.cdf_spline = UnivariateSpline(ds, self.fraction_cdf, ext=3, s=0)
                self.pdf_spline = UnivariateSpline(ds, self.fraction_cdf, ext=3, s=0).derivative(1)

            # The pdf basis integral splines will be stored here
            self.basis_integrals = {}

        def _pdf(self, d):
            return max(0.0, float(self.pdf_spline(d)))

        def _cdf(self, d):
            if d > self.d_excessive:
                # Handle spline values past 1 that decrease to zero
                return 1.0
            return max(0.0, float(self.cdf_spline(d)))

        def _pdf_basis_integral(self, d, n):
            # there are slight errors with this approach - but they are OK to
            # ignore.
            # DO NOT evaluate the first point as it leads to inf values; just set
            # it to zero
            from fluids.numerics import numpy as np
            if n not in self.basis_integrals:
                ds = np.array(self.ds[1:])
                pdf_vals = self.pdf_spline(ds)
                # n may be an integer, numpy says "Integers to negative integer powers are not allowed" if we don't make it a float
                basis_integral = ds**float(n)*pdf_vals
                if self.monotonic:
                    from scipy.interpolate import PchipInterpolator
                    self.basis_integrals[n] = PchipInterpolator(ds, basis_integral, extrapolate=True).antiderivative(1)
                else:
                    from scipy.interpolate import UnivariateSpline
                    self.basis_integrals[n] = UnivariateSpline(ds, basis_integral, ext=3, s=0).antiderivative(n=1)
            return max(float(self.basis_integrals[n](d)), 0.0)