From 005309237a9eaf2775843d4f76c3da0fc5f0d518 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Mon, 9 May 2022 14:42:53 +0200
Subject: [PATCH 01/10] feat: cli version print argument

---
 src/pygaps/cli/cli.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/pygaps/cli/cli.py b/src/pygaps/cli/cli.py
index 843fca63..12886bc5 100644
--- a/src/pygaps/cli/cli.py
+++ b/src/pygaps/cli/cli.py
@@ -71,10 +71,20 @@ def main():
         action='store_true',
         help='increase verbosity',
     )
+    prs.add_argument(
+        '--version',
+        action="store_true",
+        help="Print current version.",
+    )
 
     # Execute the parse_args() method
     args = prs.parse_args()
 
+    if args.version:
+        from importlib.metadata import version
+        print(version("pygaps"))
+        return
+
     # Read the isotherm
     if not args.iso.exists():
         raise FileNotFoundError("Path does not exist.")

From dcea649c892575df339417f2f3f77b469a9354b1 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Tue, 17 May 2022 14:13:27 +0200
Subject: [PATCH 02/10] set mpl version

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index d7f786df..15729f87 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -53,7 +53,7 @@ install_requires =
 	numpy >= 1.16.5
 	scipy >= 1.4
 	pandas
-	matplotlib
+	matplotlib >= 3.4
 	coolprop >= 6.0
 	xlrd >= 1.1
 	xlwt >= 1.3

From 4106ce1603d30b0b07a31ed97eb011d4054b1a8c Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Thu, 19 May 2022 15:47:42 +0200
Subject: [PATCH 03/10] feat: remove stale conda recipe

---
 pygaps/meta.yaml | 54 ------------------------------------------------
 1 file changed, 54 deletions(-)
 delete mode 100644 pygaps/meta.yaml

diff --git a/pygaps/meta.yaml b/pygaps/meta.yaml
deleted file mode 100644
index 1956460d..00000000
--- a/pygaps/meta.yaml
+++ /dev/null
@@ -1,54 +0,0 @@
-{% set name = "pygaps" %}
-{% set version = "4.0.1" %}
-
-package:
-  name: {{ name|lower }}
-  version: {{ version }}
-
-build:
-  skip: true  # [py<37]
-  entry_points:
-    - pygaps = pygaps.cli:main
-  script: {{ PYTHON }} -m pip install . -vv
-  number: 0
-
-source:
-  url: https://pypi.io/packages/source/{{ name[0] }}/{{ name }}/pygaps-{{ version }}.tar.gz
-  sha256: 5bdf02b497a955b19114bff9e42bce39fb133324da10a2d1c3e978b2ad452fc4
-
-requirements:
-  host:
-    - pip
-    - python
-  run:
-    - coolprop >=6.0
-    - gemmi
-    - importlib-resources  # [py<39]
-    - matplotlib-base
-    - numpy >=1.16.5
-    - openpyxl
-    - pandas
-    - python
-    - requests
-    - scipy >=1.6
-    - xlrd >=1.1
-    - xlwt >=1.3
-
-test:
-  imports:
-    - pygaps
-  commands:
-    - pip check
-    - pygaps --help
-  requires:
-    - pip
-
-about:
-  home: https://github.com/pauliacomi/pygaps
-  summary: A framework for processing adsorption data for porous materials.
-  license: MIT
-  license_file: LICENSE
-
-extra:
-  recipe-maintainers:
-    - pauliacomi

From 3f9a635615c7696a14d557f9ef4a0d935b231b04 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Mon, 27 Jun 2022 18:11:06 +0200
Subject: [PATCH 04/10] feat: error check and psd_meso improvements

---
 src/pygaps/characterisation/models_kelvin.py |  26 +-
 src/pygaps/characterisation/psd_meso.py      | 473 +++++++++++--------
 src/pygaps/graphing/calc_graphs.py           |   4 +-
 tests/characterisation/test_psd_meso.py      |  12 +-
 4 files changed, 293 insertions(+), 222 deletions(-)

diff --git a/src/pygaps/characterisation/models_kelvin.py b/src/pygaps/characterisation/models_kelvin.py
index f280058f..c09f3978 100644
--- a/src/pygaps/characterisation/models_kelvin.py
+++ b/src/pygaps/characterisation/models_kelvin.py
@@ -20,7 +20,7 @@ def get_meniscus_geometry(branch: str, pore_geometry: str):
     ----------
     branch : {'ads', 'des'}
         Branch of the isotherm used.
-    geometry : {'slit', 'cylinder', 'cylinder'}
+    geometry : {'slit', 'cylinder', 'halfopen-cylinder', 'sphere'}
         Geometry of the pore.
 
     Returns
@@ -30,23 +30,31 @@ def get_meniscus_geometry(branch: str, pore_geometry: str):
 
     """
     if branch == 'ads':
-        if pore_geometry == 'cylinder':
+        if pore_geometry == 'slit':
+            m_geometry = 'hemicylindrical'
+        elif pore_geometry == 'cylinder':
             m_geometry = 'cylindrical'
+        elif pore_geometry == 'halfopen-cylinder':
+            m_geometry = 'hemispherical'
         elif pore_geometry == 'sphere':
             m_geometry = 'hemispherical'
-        elif pore_geometry == 'slit':
-            m_geometry = 'hemicylindrical'
         else:
-            raise ParameterError("Pore geometry must be either 'cylinder', 'sphere' or 'slit'.")
+            raise ParameterError(
+                "Pore geometry must be either 'slit', 'cylinder', 'halfopen-cylinder', or 'sphere'"
+            )
     elif branch == 'des':
-        if pore_geometry == 'cylinder':
+        if pore_geometry == 'slit':
+            m_geometry = 'hemicylindrical'
+        elif pore_geometry == 'cylinder':
+            m_geometry = 'hemispherical'
+        elif pore_geometry == 'halfopen-cylinder':
             m_geometry = 'hemispherical'
         elif pore_geometry == 'sphere':
             m_geometry = 'hemispherical'
-        elif pore_geometry == 'slit':
-            m_geometry = 'hemicylindrical'
         else:
-            raise ParameterError("Pore geometry must be either 'cylinder', 'sphere' or 'slit'.")
+            raise ParameterError(
+                "Pore geometry must be either 'slit', 'cylinder', 'halfopen-cylinder', or 'sphere'"
+            )
     else:
         raise ParameterError("Adsorption branch must be either 'ads' or 'des'.")
 
diff --git a/src/pygaps/characterisation/psd_meso.py b/src/pygaps/characterisation/psd_meso.py
index 049675e3..aa1fbda1 100644
--- a/src/pygaps/characterisation/psd_meso.py
+++ b/src/pygaps/characterisation/psd_meso.py
@@ -15,15 +15,18 @@
 from pygaps.utilities.exceptions import CalculationError
 from pygaps.utilities.exceptions import ParameterError
 from pygaps.utilities.exceptions import pgError
+from pygaps import logger
 
 _MESO_PSD_MODELS = ['pygaps-DH', 'BJH', 'DH']
-_PORE_GEOMETRIES = ['slit', 'cylinder', 'sphere']
+_PORE_GEOMETRIES = ['slit', 'cylinder', 'halfopen-cylinder', 'sphere']
+_MENISCUS_GEOMETRIES = ["hemicylindrical", "cylindrical", "hemispherical"]
 
 
 def psd_mesoporous(
     isotherm: "PointIsotherm | ModelIsotherm",
     psd_model: str = 'pygaps-DH',
     pore_geometry: str = 'cylinder',
+    meniscus_geometry: str = None,
     branch: str = 'des',
     thickness_model:
     "str | PointIsotherm | ModelIsotherm | t.Callable[[float], float]" = 'Harkins/Jura',
@@ -35,6 +38,9 @@ def psd_mesoporous(
     Calculate the mesopore size distribution using a Kelvin-type model.
 
     Expected pore geometry must be specified as ``pore_geometry``.
+    If the meniscus (adsorbed/gas interface) geometry is known, it
+    can be equally specified, otherwise it will be inferred from the
+    pore geometry.
 
     The ``thickness_model`` parameter is a string which names the thickness
     equation which should be used. Alternatively, a user can implement their own
@@ -49,6 +55,8 @@ def psd_mesoporous(
     psd_model : str
         The pore size distribution model to use.
     pore_geometry : str
+        The geometry of the meniscus (adsorbed/gas interface) in the pores.
+    meniscus_geometry : str, optional
         The geometry of the adsorbent pores.
     branch : {'ads', 'des'}, optional
         Branch of the isotherm to use. It defaults to desorption.
@@ -58,18 +66,23 @@ def psd_mesoporous(
         A custom user kelvin model. It should be a callable that only takes
         relative pressure as an argument.
     p_limits : tuple[float, float]
-        Pressure range in which to calculate PSD, defaults to entire isotherm.
+        Pressure range in which to calculate PSD, defaults to (0.1, 0.99).
     verbose : bool
         Prints out extra information on the calculation and graphs the results.
 
     Returns
     -------
     dict
-        A dictionary with the pore widths and the pore distributions, of the form:
+        A dictionary of results of the form:
 
-        - ``pore_widths`` (array) : the widths of the pores
-        - ``pore_distribution`` (array) : contribution of each pore width to the
-          overall pore distribution
+        - ``pore_widths`` (ndarray) : the widths (or diameter) of the pores, nm
+        - ``pore_volumes`` (ndarray) : pore volume for each pore width, cm3/material
+        - ``pore_volume_cumulative`` (ndarray) : cumulative sum of pore volumes, in cm3/material
+        - ``pore_distribution`` (ndarray) : contribution of each pore width to the
+          overall pore distribution, cm3/material/nm
+        - ``pore_areas`` (ndarray) : specific area for each pore width, m2/material
+        - ``pore_area_total`` (float) : total specific area, m2/material
+        - ``limits`` (tuple[int, int]) : indices of selection limits
 
     Raises
     ------
@@ -82,36 +95,47 @@ def psd_mesoporous(
     -----
     Calculates the pore size distribution using a 'classical' model which attempts to
     describe the adsorption in a pore as a combination of a statistical thickness and
-    a condensation/evaporation behaviour described by surface tension.
+    a condensation/evaporation behaviour related to adsorbate surface tension.
     It is based on solving the following equation:
 
     .. math::
 
-        \Delta V_n = R V_p + \Delta t \sum_{i=1}^{n-1} A_i
+        \Delta V_n = V_{k,n} + V_{t,n}
 
-    Which states that the volume adsorbed or desorbed during an isotherm step
-    is a sum of the capillary condensation / evaporation volume, expressed as a
-    function of the total pore volume, and the volume of the thickness increase
-    or decrease in a pore.
+    Which states that the volume adsorbed or desorbed during a pressure step,
+    :math:`\Delta V_n`, can be described as a sum of the volume involved in
+    capillary condensation / evaporation (:math:`\Delta V_{k,n}`), and the
+    volume corresponding to the increase / decrease of the adsorbed layer
+    thickness in the surface of all non-filled pores (:math:`\Delta V_{t,n}`).
+
+    Expressions are derived for the pore volume as a function of pore geometry,
+    the shape of the liquid-gas interface (meniscus), relationship between
+    pore width and critical condensation width (:math:`R_{p}`), rearranging
+    the equation to:
+
+    .. math::
+
+        V_{p,n} = (\Delta V_n - V_{t,n}) R_p
 
     Currently, the methods provided are:
 
-    - the pygaps-DH model, an expanded DH model for multiple geometries
+    - the pygaps-DH model, a custom expanded DH model for multiple pore geometries
     - the original BJH or Barrett, Joyner and Halenda method
-    - the original DH or Dollimore-Heal method, an extension of the BJH method
+    - the original DH or Dollimore-Heal method
 
     According to Rouquerol [#]_, in adopting this approach, it is assumed that:
 
     - The Kelvin equation is applicable over the pore range (mesopores).
-      Therefore in pores which are below a certain size (around 2.5 nm), the
-      granularity of the liquid-vapour interface becomes too large for classical
-      bulk methods to be applied.
+      In pores which are below a certain size (around 2.5 nm), the
+      granularity of the liquid-vapour interface becomes too large the method
+      to hold.
     - The meniscus curvature is controlled be the pore size and shape. Ideal
       shapes for the curvature are assumed.
     - The pores are rigid and of well defined shape. They are considered
-      open-ended and non-intersecting
-    - The filling/emptying of each pore does not depend on its location.
-    - The adsorption on the pore walls is not different from surface adsorption.
+      open-ended and non-intersecting.
+    - The filling/emptying of each pore does not depend on its location,
+      i.e. no diffusion or blocking effects.
+    - Adsorption on the pore walls is not different from surface adsorption.
 
     .. caution::
 
@@ -142,9 +166,14 @@ def psd_mesoporous(
         )
     if pore_geometry not in _PORE_GEOMETRIES:
         raise ParameterError(
-            f"Geometry {pore_geometry} not an option for pore size distribution. ",
+            f"Geometry {pore_geometry} not an option for pore size distribution. "
             f"Available geometries are {_PORE_GEOMETRIES}"
         )
+    if meniscus_geometry and meniscus_geometry not in _MENISCUS_GEOMETRIES:
+        raise ParameterError(
+            f"Meniscus geometry {meniscus_geometry} not an option for pore size distribution. "
+            f"Available geometries are {_MENISCUS_GEOMETRIES}"
+        )
     if branch not in ['ads', 'des']:
         raise ParameterError(
             f"Branch '{branch}' not an option for PSD.", "Select either 'ads' or 'des'"
@@ -158,12 +187,12 @@ def psd_mesoporous(
     surface_tension = isotherm.adsorbate.surface_tension(isotherm.temperature)
 
     # Read data in, depending on branch requested
-    loading = isotherm.loading(
+    volume_adsorbed = isotherm.loading(
         branch=branch,
-        loading_basis='molar',
-        loading_unit='mmol',
+        loading_basis='volume_liquid',
+        loading_unit='cm3',
     )
-    if loading is None:
+    if volume_adsorbed is None:
         raise ParameterError("The isotherm does not have the required branch for this calculation")
     try:
         pressure = isotherm.pressure(
@@ -178,8 +207,8 @@ def psd_mesoporous(
 
     # If on an desorption branch, data will be reversed
     if branch == 'des':
-        loading = loading[::-1]
         pressure = pressure[::-1]
+        volume_adsorbed = volume_adsorbed[::-1]
 
     # select the maximum and minimum of the points and the pressure associated
     minimum = 0
@@ -187,7 +216,7 @@ def psd_mesoporous(
 
     # Set default values
     if p_limits is None:
-        p_limits = (None, 0.95)
+        p_limits = (0.1, 0.99)
 
     if p_limits[0]:
         minimum = numpy.searchsorted(pressure, p_limits[0])
@@ -199,55 +228,59 @@ def psd_mesoporous(
             "in the selected region."
         )
     pressure = pressure[minimum:maximum + 1]
-    loading = loading[minimum:maximum + 1]
+    volume_adsorbed = volume_adsorbed[minimum:maximum + 1]
 
-    # calculated volume adsorbed
-    volume_adsorbed = loading * molar_mass / liquid_density / 1000
-
-    # Thickness model definitions
+    # Thickness model
     t_model = get_thickness_model(thickness_model)
 
-    # Kelvin model definitions
-    k_model_args = {
-        "meniscus_geometry": get_meniscus_geometry(branch, pore_geometry),
-        "temperature": isotherm.temperature,
-        "liquid_density": liquid_density,
-        "adsorbate_molar_mass": molar_mass,
-        "adsorbate_surface_tension": surface_tension
-    }
-    k_model = get_kelvin_model(kelvin_model, **k_model_args)
+    # Kelvin model
+    if not meniscus_geometry:
+        meniscus_geometry = get_meniscus_geometry(branch, pore_geometry)
+    k_model = get_kelvin_model(
+        kelvin_model,
+        meniscus_geometry=meniscus_geometry,
+        temperature=isotherm.temperature,
+        liquid_density=liquid_density,
+        adsorbate_molar_mass=molar_mass,
+        adsorbate_surface_tension=surface_tension
+    )
 
     # Call specified pore size distribution function
     if psd_model == 'pygaps-DH':
-        pore_widths, pore_dist, pore_vol_cum = psd_pygapsdh(
-            volume_adsorbed, pressure, pore_geometry, t_model, k_model
-        )
+        results = psd_pygapsdh(volume_adsorbed, pressure, pore_geometry, t_model, k_model)
     elif psd_model == 'BJH':
-        pore_widths, pore_dist, pore_vol_cum = psd_bjh(
-            volume_adsorbed, pressure, pore_geometry, t_model, k_model
-        )
+        results = psd_bjh(volume_adsorbed, pressure, pore_geometry, t_model, k_model)
     elif psd_model == 'DH':
-        pore_widths, pore_dist, pore_vol_cum = psd_dollimore_heal(
-            volume_adsorbed, pressure, pore_geometry, t_model, k_model
+        results = psd_dollimore_heal(volume_adsorbed, pressure, pore_geometry, t_model, k_model)
+
+    results['pore_volume_cumulative'] = numpy.cumsum(results['pore_volumes'])
+    results['pore_volume_cumulative'] = (
+        results['pore_volume_cumulative'] - results['pore_volume_cumulative'][-1] +
+        volume_adsorbed[-1]
+    )
+    results['limits'] = (minimum, maximum)
+    results['pore_area_total'] = sum(results['pore_areas'])
+
+    if any(results['pore_volume_cumulative'] < 0):
+        logger.warning(
+            "Negative values encountered in cumulative pore volumes. "
+            "It is very likely that the model or its limits are wrong. "
+            "Check that your pore geometry, meniscus geometry and thickness function "
+            "are suitable for the material."
         )
 
     # Plot if verbose
     if verbose:
         from pygaps.graphing.calc_graphs import psd_plot
         psd_plot(
-            pore_widths,
-            pore_dist,
-            pore_vol_cum=pore_vol_cum,
+            results["pore_widths"],
+            results["pore_distribution"],
+            pore_vol_cum=results['pore_volume_cumulative'],
             method=psd_model,
             left=1.5,
         )
 
-    return {
-        'pore_widths': pore_widths,
-        'pore_distribution': pore_dist,
-        'pore_volume_cumulative': pore_vol_cum,
-        'limits': (minimum, maximum),
-    }
+    return results
 
 
 def psd_pygapsdh(
@@ -265,25 +298,27 @@ def psd_pygapsdh(
     Parameters
     ----------
     volume_adsorbed : array
-        Volume adsorbed of "liquid" phase in cm3.
+        Volume adsorbed of "liquid" phase in cm3/material.
     relative_pressure : array
         Relative pressure.
     pore_geometry : str
         The geometry of the pore, eg. 'sphere', 'cylinder' or 'slit'.
     thickness_model : callable
         Function which returns the thickness of the adsorbed layer
-        at a pressure p.
+        at a pressure p, in nm.
     condensation_model : callable
-        Function which returns the critical kelvin radius at a pressure p.
+        Function which returns the critical kelvin radius at
+        a pressure p, in nm.
 
     Returns
     -------
-    pore widths : array
-        Widths of the pores.
-    pore_dist : array
-        Amount of each pore width.
-    pore_vol_cum : array
-        Cumulative pore volume.
+    dict
+        A dictionary of results of the form:
+
+        - ``pore_widths`` (ndarray) : the widths (or diameter) of the pores, nm
+        - ``pore_volumes`` (ndarray) : pore volume for each pore width, cm3/material
+        - ``pore_areas`` (ndarray) : specific area for each pore width, m2/material
+        - ``pore_distribution`` (ndarray) : volumetric pore distribution, cm3/material/nm
 
     Notes
     -----
@@ -297,7 +332,7 @@ def psd_pygapsdh(
 
     .. math::
 
-        w_{p} = 2 r_k + 2 t
+        w_p = 2 * r_p = (r_k + t)
 
     The Kelvin radius and layer thickness are calculated based on the
     functions provided, no assessment of the meniscus shape or thickness is
@@ -307,9 +342,9 @@ def psd_pygapsdh(
 
     .. math::
 
-        V_p = \Big(\Delta V_n - \Delta t_n \sum_{i=1}^{n-1}
-                \Big(\frac{\bar{w}_{p,i} - 2 t_{n,i}}{\bar{w}_{p,i}}\Big)^{(l-1)}
-                \frac{2 l \Delta V_{p,i}}{w_{p,i}}\Big)
+        V_p = \Big[\Delta V_n - \Delta t_n \sum_{i=1}^{n-1}
+                \Big(\frac{\bar{w}_{p,i} - 2 t_{n}}{\bar{w}_{p,i}}\Big)^{l-1}
+                \frac{2 l V_{p,i}}{w_{p,i}}\Big]
                 \Big(\frac{\bar{w}_p}{\bar{w}_p - 2 t_n}\Big)^l
 
     Where :
@@ -318,11 +353,11 @@ def psd_pygapsdh(
         - :math:`\bar{r}_p` is the average pore radius calculated as a sum of the
           kelvin radius and layer thickness of the pores at pressure p between two
           measurement points
-        - :math:`\bar{r}_k` is the average kelvin radius between two measurement points
+        - :math:`\bar{r}_k` is the average kelvin radius between two pressure points
         - :math:`t_n` is the layer thickness at point n
-        - :math:`\bar{t}_n` is the average layer thickness between two measurement points
-        - :math:`\Delta t_n` is the change in layer thickness between two measurement points
-        - :math:`l` is a characteristic dimension of the system
+        - :math:`\bar{t}_n` is the average layer thickness between two pressure points
+        - :math:`\Delta t_n` is the change in layer thickness between two pressure points
+        - :math:`l` is the characteristic dimension of the system
 
     In order to account for different geometries, factors are calculated in
     terms of a characteristic number of the system:
@@ -362,6 +397,7 @@ def psd_pygapsdh(
         raise ParameterError("The length of the pressure and loading arrays do not match.")
 
     # Pore geometry specifics
+    # A constant is defined which is used in the general formula.
     if pore_geometry == 'slit':
         c_length = 1
     elif pore_geometry == 'cylinder':
@@ -372,51 +408,57 @@ def psd_pygapsdh(
         raise ParameterError("Unknown pore geometry.")
 
     # We reverse the arrays, starting from the highest loading
-    volume_adsorbed = volume_adsorbed[::-1]
-    relative_pressure = relative_pressure[::-1]
+    volume_adsorbed = volume_adsorbed[::-1]  # [cm3/mat]
+    relative_pressure = relative_pressure[::-1]  # [unitless]
 
-    # Calculate the adsorbed volume of liquid and diff
+    # Calculate the adsorbed volume of liquid and diff, [cm3/mat]
     d_volume = -numpy.diff(volume_adsorbed)
 
-    # Generate the thickness curve, average and diff
+    # Generate the thickness curve, average and diff, [nm]
     thickness = thickness_model(relative_pressure)
     avg_thickness = (thickness[:-1] + thickness[1:]) / 2
     d_thickness = -numpy.diff(thickness)
 
-    # Generate the Kelvin pore radii and average
-    kelvin_radius = condensation_model(relative_pressure)
+    # Generate the Kelvin pore radii and average, [nm]
+    kelvin_radii = condensation_model(relative_pressure)
 
     # Critical pore radii as a combination of the adsorbed
-    # layer thickness and kelvin pore radius, with average and diff
-    pore_widths = 2 * (thickness + kelvin_radius)
+    # layer thickness and kelvin pore radius, with average and diff, [nm]
+    pore_widths = 2 * (thickness + kelvin_radii)
     avg_pore_widths = (pore_widths[:-1] + pore_widths[1:]) / 2
     d_pore_widths = -numpy.diff(pore_widths)
 
+    # Calculate the ratios of the pore to the evaporated capillary "core", [unitless]
+    ratio_factors = (avg_pore_widths / (avg_pore_widths - 2 * avg_thickness))**2
+
     # Now we can iteratively calculate the pore size distribution
-    sum_area_factor = 0
-    pore_volumes = []
+    sum_area_correction = 0  # cm3/nm
+    pore_areas = numpy.zeros_like(avg_pore_widths)  # areas of pore populations [m2/mat]
+    pore_volumes = numpy.zeros_like(avg_pore_widths)  # volume of pore populations, [cm3/mat]
 
     for i, avg_pore_width in enumerate(avg_pore_widths):
 
-        # Calculate the ratio of the pore to the evaporated capillary "core"
-        ratio_factor = (avg_pore_width / (avg_pore_width - 2 * thickness[i]))**c_length
-
-        # Calculate the volume desorbed from thinning of all pores previously emptied
-        thickness_factor = -d_thickness[i] * sum_area_factor
+        # Volume desorbed from thinning of all pores previously emptied, [cm3]
+        d_thickness_volume = d_thickness[i] * sum_area_correction
 
-        # Equation for pore volume, then add to the array
-        pore_volume = (d_volume[i] + thickness_factor) * ratio_factor
-        pore_volumes.append(pore_volume)
+        # Volume of newly emptied pore, [cm3]
+        pore_volume = (d_volume[i] - d_thickness_volume) * ratio_factors[i]
+        pore_volumes[i] = pore_volume
 
-        # Calculate the area of the newly emptied pore then add it to the total pore area
-        pore_area_correction = ((avg_pore_width - 2 * avg_thickness[i]) /
-                                avg_pore_width)**(c_length - 1)
-        pore_avg_area = 2 * c_length * pore_volume / avg_pore_width
-        sum_area_factor += pore_area_correction * pore_avg_area
+        # Area of the newly emptied pore, [m2]
+        pore_geometry_correction = ((avg_pore_width - 2 * avg_thickness[i]) /
+                                    avg_pore_width)**(c_length - 1)
+        pore_area = 2 * c_length * pore_volume / avg_pore_width  # cm3/nm
+        pore_areas[i] = pore_area * 1e3  # cm3/nm = 1e-6 m3/ 1e-9m
 
-    pore_dist = pore_volumes / d_pore_widths
+        sum_area_correction += pore_geometry_correction * pore_area
 
-    return pore_widths[:0:-1], pore_dist[::-1], numpy.cumsum(pore_volumes[::-1])
+    return {
+        "pore_widths": pore_widths[:0:-1],  # [nm]
+        "pore_areas": pore_areas[::-1],  # [m2/mat]
+        "pore_volumes": pore_volumes[::-1],  # [cm3/mat]
+        "pore_distribution": (pore_volumes / d_pore_widths)[::-1],  # [cm3/mat/nm]
+    }
 
 
 def psd_bjh(
@@ -427,48 +469,51 @@ def psd_bjh(
     condensation_model: t.Callable,
 ):
     r"""
-    Calculate a pore size distribution using the BJH method.
+    Calculate a pore size distribution using the original BJH method.
     This function should not be used with isotherms (use instead
     :func:`pygaps.characterisation.psd_meso.psd_mesoporous`).
 
     Parameters
     ----------
     volume_adsorbed : array
-        Volume adsorbed of "liquid" phase in cm3.
+        Volume adsorbed of "liquid" phase in cm3/material.
     relative_pressure : array
         Relative pressure.
     pore_geometry : str
         The geometry of the pore, eg. 'sphere', 'cylinder' or 'slit'.
     thickness_model : callable
         Function which returns the thickness of the adsorbed layer
-        at a pressure p.
+        at a pressure p, in nm.
     condensation_model : callable
-        Function which returns the critical kelvin radius at a pressure p.
+        Function which returns the critical kelvin radius at
+        a pressure p, in nm.
 
     Returns
     -------
-    pore widths : array
-        Widths of the pores.
-    pore_dist : array
-        Amount of each pore width.
-    pore_vol_cum : array
-        Cumulative pore volume.
+    dict
+        A dictionary of results of the form:
+
+        - ``pore_widths`` (ndarray) : the widths (or diameter) of the pores, nm
+        - ``pore_volumes`` (ndarray) : pore volume for each pore width, cm3/material
+        - ``pore_areas`` (ndarray) : specific area for each pore width, m2/material
+        - ``pore_distribution`` (ndarray) : volumetric pore distribution, cm3/material/nm
 
     Notes
     -----
+
     The BJH or Barrett, Joyner and Halenda [#]_ method for calculation of pore
     size distribution is based on a classical description of the adsorbate
     behaviour in the adsorbent pores. Under this method, the adsorbate is
-    adsorbing on the pore walls in a predictable way, and decreasing the
-    apparent pore volume until condensation takes place, filling the entire
-    pore. The two variables, layer thickness and radius where condensation takes
-    place can be modelled by a thickness model (such as Halsey, Harkins & Jura,
-    etc.) and a critical radius model for condensation/evaporation, based on a
-    form of the Kelvin equation.
+    adsorbing on the pore walls forming a layer of known thickness, therefore
+    decreasing the apparent pore volume until condensation takes place, filling
+    the entire pore. The two variables, layer thickness and critical pore width
+    where condensation takes place can be respectively modelled by a thickness
+    model (such as Halsey, Harkins & Jura, etc.) and a model for
+    condensation/evaporation based on a form of the Kelvin equation.
 
     .. math::
 
-        r_p = t + r_k
+        1/2 w_p = r_p = r_k + t
 
     The original model used the desorption curve as a basis for calculating pore
     size distribution. Between two points of the curve, the volume desorbed can
@@ -479,8 +524,8 @@ def psd_bjh(
 
     .. math::
 
-        V_p = \Big(\Delta V_n - \Delta t_n \sum_{i=1}^{n-1}
-                \Big(\frac{\bar{r}_{p,i} - t_{n,i}}{\bar{r}_{p,i}}\Big) \frac{2 \Delta V_{p,i}}{r_{p,i}}\Big)
+        V_p = \Big[\Delta V_n - \Delta t_n \sum_{i=1}^{n-1}
+                \Big(\frac{\bar{r}_{p,i} - t_{n}}{\bar{r}_{p,i}}\Big) A_{p,i}\Big]
                 \Big(\frac{\bar{r}_p}{\bar{r}_k + \Delta t_n}\Big)^2
 
     Where :
@@ -515,7 +560,8 @@ def psd_bjh(
     if len(volume_adsorbed) != len(relative_pressure):
         raise ParameterError("The length of the pressure and loading arrays do not match.")
 
-    if pore_geometry in ('slit', 'sphere'):
+    # Pore geometry specifics
+    if pore_geometry != 'cylinder':
         raise ParameterError(
             "The BJH method is provided for compatibility and only applicable"
             " to cylindrical pores. Use the pyGAPS-DH method for other options."
@@ -525,48 +571,57 @@ def psd_bjh(
     volume_adsorbed = volume_adsorbed[::-1]
     relative_pressure = relative_pressure[::-1]
 
-    # Calculate the adsorbed volume of liquid and diff
+    # Calculate the adsorbed volume of liquid and diff, [cm3/mat]
     d_volume = -numpy.diff(volume_adsorbed)
 
-    # Generate the thickness curve, average and diff
+    # Generate the thickness curve, average and diff, [nm]
     thickness = thickness_model(relative_pressure)
     avg_thickness = (thickness[:-1] + thickness[1:]) / 2
     d_thickness = -numpy.diff(thickness)
 
-    # Generate the Kelvin pore radii and average
-    kelvin_radius = condensation_model(relative_pressure)
-    avg_k_radius = (kelvin_radius[:-1] + kelvin_radius[1:]) / 2
+    # Generate the Kelvin pore radii and average, [nm]
+    kelvin_radii = condensation_model(relative_pressure)
+    avg_k_radii = (kelvin_radii[:-1] + kelvin_radii[1:]) / 2
 
     # Critical pore radii as a combination of the adsorbed
-    # layer thickness and kelvin pore radius, with average and diff
-    pore_widths = 2 * (thickness + kelvin_radius)
-    avg_pore_widths = 2 * (avg_thickness + avg_k_radius)
-    d_pore_widths = -numpy.diff(pore_widths)
+    # layer thickness and kelvin pore radius, with average and diff, [nm]
+    pore_radii = thickness + kelvin_radii
+    avg_pore_radii = (pore_radii[:-1] + pore_radii[1:]) / 2
+    d_pore_radii = -numpy.diff(pore_radii)
 
-    # Now we can iteratively calculate the pore size distribution
-    sum_area_factor = 0
-    pore_volumes = []
+    # Calculate the ratio of the pore to the evaporated capillary "core", [unitless]
+    ratio_factors = (avg_pore_radii / (avg_k_radii + d_thickness))**2
 
-    for i, avg_pore_width in enumerate(avg_pore_widths):
+    # Now we can iteratively calculate the pore size distribution
+    pore_areas = numpy.zeros_like(avg_pore_radii)  # areas of pore populations [m2/mat]
+    pore_volumes = numpy.zeros_like(avg_pore_radii)  # volume of pore populations, [cm3/mat]
 
-        # Calculate the ratio of the pore to the evaporated capillary "core"
-        ratio_factor = (avg_pore_width / (2 * (avg_k_radius[i] + d_thickness[i])))**2
+    for i, avg_pore_rad in enumerate(avg_pore_radii):
 
-        # Calculate the volume desorbed from thinning of all pores previously emptied
-        thickness_factor = -d_thickness[i] * sum_area_factor
+        # Pore area correction
+        sum_area_factor = 0  # [m2]
+        for x in range(i):
+            area_ratio = (avg_pore_radii[x] - avg_thickness[i]) / avg_pore_radii[x]
+            sum_area_factor += area_ratio * pore_areas[x]
 
-        # Equation for pore volume, then add to the array
-        pore_volume = (d_volume[i] + thickness_factor) * ratio_factor
-        pore_volumes.append(pore_volume)
+        # Calculate the volume desorbed from thinning of all pores previously emptied, [cm3/mat]
+        # nm * m3 = 1e-7 cm * 1e4 cm2 = 1e-3 cm3
+        d_thickness_volume = d_thickness[i] * sum_area_factor * 1e-3
 
-        # Calculate the area of the newly emptied pore then add it to the total pore area
-        pore_area_correction = (avg_pore_width - 2 * avg_thickness[i]) / avg_pore_width
-        pore_avg_area = (4 * pore_volume / avg_pore_width)
-        sum_area_factor += pore_area_correction * pore_avg_area
+        # Calculate volume of newly emptied pore, [cm3]
+        pore_volume = (d_volume[i] - d_thickness_volume) * ratio_factors[i]
+        pore_volumes[i] = pore_volume
 
-    pore_dist = pore_volumes / d_pore_widths
+        # Calculate the area of the newly emptied pore, [m2]
+        pore_area = 2 * pore_volume / avg_pore_rad * 1e3  # cm3/nm = 1e-6 m3/ 1e-9m
+        pore_areas[i] = pore_area
 
-    return pore_widths[:0:-1], pore_dist[::-1], numpy.cumsum(pore_volumes[::-1])
+    return {
+        "pore_widths": pore_radii[:0:-1] * 2,  # [nm]
+        "pore_areas": pore_areas[::-1],  # [m2/mat]
+        "pore_volumes": pore_volumes[::-1],  # [cm3/mat]
+        "pore_distribution": (pore_volumes / d_pore_radii / 2)[::-1],  # [cm3/mat/nm]
+    }
 
 
 def psd_dollimore_heal(
@@ -577,51 +632,52 @@ def psd_dollimore_heal(
     condensation_model: t.Callable,
 ):
     r"""
-    Calculate a pore size distribution using the Dollimore-Heal method.
+    Calculate a pore size distribution using the original Dollimore-Heal method.
     This function should not be used with isotherms (use instead
     :func:`pygaps.characterisation.psd_meso.psd_mesoporous`).
 
     Parameters
     ----------
     volume_adsorbed : array
-        Volume adsorbed of "liquid" phase in cm3.
+        Volume adsorbed of "liquid" phase in cm3/material.
     relative_pressure : array
         Relative pressure.
     pore_geometry : str
         The geometry of the pore, eg. 'sphere', 'cylinder' or 'slit'.
     thickness_model : callable
         Function which returns the thickness of the adsorbed layer
-        at a pressure p.
+        at a pressure p, in nm.
     condensation_model : callable
-        Function which returns the critical kelvin radius at a pressure p.
+        Function which returns the critical kelvin radius at
+        a pressure p, in nm.
 
     Returns
     -------
-    pore widths : array
-        Widths of the pores.
-    pore_dist : array
-        Amount of each pore width.
-    pore_vol_cum : array
-        Cumulative pore volume.
+    dict
+        A dictionary of results of the form:
+
+        - ``pore_widths`` (ndarray) : the widths (or diameter) of the pores, nm
+        - ``pore_volumes`` (ndarray) : pore volume for each pore width, cm3/material
+        - ``pore_areas`` (ndarray) : specific area for each pore width, m2/material
+        - ``pore_distribution`` (ndarray) : volumetric pore distribution, cm3/material/nm
 
     Notes
     -----
     The DH or Dollimore-Heal method [#]_ of calculation of pore size distribution is an
-    extension of the BJH method, which takes into account the geometry of the pores
-    by introducing a length component.
+    extension of the BJH method.
 
     Like the BJH method, it is based on a classical description of the adsorbate
     behaviour in the adsorbent pores. Under this method, the adsorbate is
     adsorbing on the pore walls in a predictable way, and decreasing the
-    apparent pore volume until condensation takes place, filling the entire
-    pore. The two variables, layer thickness and radius where condensation takes
-    place can be modelled by a thickness model (such as Halsey, Harkins & Jura,
-    etc.) and a critical radius model for condensation/evaporation, based on a
-    form of the Kelvin equation.
+    apparent pore radius until condensation takes place, filling the entire
+    pore. The two components, layer thickness (t) and radius (r_k) where
+    condensation takes place can be modelled by a thickness model (such as
+    Halsey, Harkins & Jura, etc.) and a critical radius model for
+    condensation/evaporation, based on a form of the Kelvin equation.
 
     .. math::
 
-        r_p = t + r_k
+        1/2 w_p = r_p = r_k + t
 
     The original model used the desorption curve as a basis for calculating pore
     size distribution. Between two points of the curve, the volume desorbed can
@@ -632,17 +688,18 @@ def psd_dollimore_heal(
 
     .. math::
 
-        V_p = \Big(\Delta V_n - \Delta t_n \sum_{i=1}^{n-1} \frac{2 \Delta V_{p,i}}{r_{p,i}}
-                + \Delta t_n \bar{t}_n \sum_{i=1}^{n-1} \frac{2 \Delta V_{p,i}}{r_{p,i}^2}
-                \Big)\Big(\frac{\bar{r}_p}{\bar{r}_p - t_n}\Big)^2
+        V_p = \Big(\Delta V_n - \Delta V_m\Big)\Big(\frac{\bar{r}_p}{\bar{r}_p - t_n}\Big)^2
+
+        V_p = \Big[\Delta V_n - \Delta t_n \sum_{i=1}^{n-1} A_{p,i}
+                + \Delta t_n \bar{t}_n \sum_{i=1}^{n-1} 2 \pi L_{p,i}
+                \Big]\Big(\frac{\bar{r}_p}{\bar{r}_p - t_n}\Big)^2
 
     Where :
-        - :math:`V_p` is the volume of pores of size :math:`\bar{r}_p`
-        - :math:`\Delta V_n` is the adsorbed volume change between two points
         - :math:`\bar{r}_p` is the average pore radius calculated as a sum of the
           kelvin radius and layer thickness of the pores at pressure p between two
           measurement points
-        - :math:`t_n` is the layer thickness at point n
+        - :math:`V_p` is the volume of pores of size :math:`\bar{r}_p`
+        - :math:`\Delta V_n` is the adsorbed volume change between two points
         - :math:`\bar{t}_n` is the average layer thickness between two measurement points
         - :math:`\Delta t_n` is the change in layer thickness between two measurement points
 
@@ -663,58 +720,68 @@ def psd_dollimore_heal(
     if len(volume_adsorbed) != len(relative_pressure):
         raise ParameterError("The length of the pressure and loading arrays do not match.")
 
-    if pore_geometry in ('slit', 'sphere'):
+    # Pore geometry specifics
+    if pore_geometry != 'cylinder':
         raise ParameterError(
-            "The DH method is provided for compatibility and only applicable"
-            " to cylindrical pores. Use the pyGAPS-DH method for other options."
+            "The DH method is provided for compatibility and only applicable "
+            "to cylindrical pores. Use the pyGAPS-DH method for other options."
         )
 
     # We reverse the arrays, starting from the highest loading
-    volume_adsorbed = volume_adsorbed[::-1]
-    relative_pressure = relative_pressure[::-1]
+    volume_adsorbed = volume_adsorbed[::-1]  # [cm3/mat]
+    relative_pressure = relative_pressure[::-1]  # [unitless]
 
-    # Calculate the first differential of volume adsorbed
+    # Calculate the first differential of volume adsorbed, [cm3/mat]
     d_volume = -numpy.diff(volume_adsorbed)
 
-    # Generate the thickness curve, average and diff
+    # Generate the thickness curve, average and diff, [nm]
     thickness = thickness_model(relative_pressure)
     avg_thickness = (thickness[:-1] + thickness[1:]) / 2
     d_thickness = -numpy.diff(thickness)
 
-    # Generate the Kelvin pore radii and average
-    kelvin_radius = condensation_model(relative_pressure)
-    avg_k_radius = (kelvin_radius[:-1] + kelvin_radius[1:]) / 2
+    # Generate the Kelvin pore radii and average, [nm]
+    kelvin_radii = condensation_model(relative_pressure)
+    avg_k_radii = (kelvin_radii[:-1] + kelvin_radii[1:]) / 2
 
     # Critical pore radii as a combination of the adsorbed
-    # layer thickness and kelvin pore radius, with average and diff
-    pore_widths = 2 * (thickness + kelvin_radius)
-    avg_pore_widths = 2 * (avg_thickness + avg_k_radius)
-    d_pore_widths = -numpy.diff(pore_widths)
+    # layer thickness and kelvin pore radius, with average and diff, [nm]
+    pore_radii = thickness + kelvin_radii
+    avg_pore_radii = (pore_radii[:-1] + pore_radii[1:]) / 2
+    d_pore_radii = -numpy.diff(pore_radii)
 
-    # Now we can iteratively calculate the pore size distribution
-    sum_area_factor = 0
-    sum_length_factor = 0
-    pore_volumes = []
+    # Calculate the ratios of the pore to the evaporated capillary "core", [unitless]
+    ratio_factors = (avg_pore_radii / (avg_k_radii + d_thickness))**2
 
-    for i, avg_pore_width in enumerate(avg_pore_widths):
+    # Now we can iteratively calculate the pore size distribution
+    sum_area_factor = 0  # cm3/nm
+    sum_2pi_length_factor = 0  # cm3/nm2
+    pore_areas = numpy.zeros_like(avg_pore_radii)  # areas of pore populations [m2/mat]
+    pore_volumes = numpy.zeros_like(avg_pore_radii)  # volume of pore populations, [cm3/mat]
 
-        # Calculate the ratio of the pore to the evaporated capillary "core"
-        ratio_factor = (avg_pore_width / (avg_pore_width - 2 * thickness[i]))**2
+    for i, avg_pore_radius in enumerate(avg_pore_radii):
 
         # Calculate the volume desorbed from thinning of all pores previously emptied
-        thickness_factor = - d_thickness[i] * sum_area_factor + \
-            d_thickness[i] * avg_thickness[i] * sum_length_factor
+        # dt * \sum A - t * dt * 2 * \pi * \sum L
+        d_thickness_volume = d_thickness[i] * sum_area_factor - \
+            d_thickness[i] * avg_thickness[i] * sum_2pi_length_factor # [cm3/mat]
 
-        # Equation for pore volume, then add to the array
-        pore_volume = (d_volume[i] + thickness_factor) * ratio_factor
-        pore_volumes.append(pore_volume)
+        # Pore volume, then store
+        # dVp = (dV - dVt) * Rp
+        pore_volume = (d_volume[i] - d_thickness_volume) * ratio_factors[i]  # [cm3/mat]
+        pore_volumes[i] = pore_volume
 
-        # Calculate the two factors in the DH method, for area and length
-        pore_avg_area = (4 * pore_volume / avg_pore_width)
-        pore_avg_length = (8 * pore_volume / avg_pore_width**2)
-        sum_area_factor += pore_avg_area
-        sum_length_factor += pore_avg_length
+        # Calculate the two correction factors in the DH method, for area and length
+        # Ap = 2 * dVp / rp
+        pore_area = 2 * pore_volume / avg_pore_radius  # cm3/nm
+        pore_areas[i] = pore_area * 1e3  # cm3/nm = 1e-6 m3/ 1e-9m = 1e3 m2
+        sum_area_factor += pore_area
 
-    pore_dist = pore_volumes / d_pore_widths
+        # 2 * \pi * Lp = Ap / rp
+        sum_2pi_length_factor += pore_area / avg_pore_radius  # cm3/nm2
 
-    return pore_widths[:0:-1], pore_dist[::-1], numpy.cumsum(pore_volumes[::-1])
+    return {
+        "pore_widths": pore_radii[:0:-1] * 2,  # [nm]
+        "pore_areas": pore_areas[::-1],  # [m2/mat]
+        "pore_volumes": pore_volumes[::-1],  # [cm3/mat]
+        "pore_distribution": (pore_volumes / d_pore_radii / 2)[::-1],  # [cm3/mat/nm]
+    }
diff --git a/src/pygaps/graphing/calc_graphs.py b/src/pygaps/graphing/calc_graphs.py
index cdbe7e8e..e2ec8af5 100644
--- a/src/pygaps/graphing/calc_graphs.py
+++ b/src/pygaps/graphing/calc_graphs.py
@@ -428,9 +428,9 @@ def formatter(x, pos):
 
     ax.set_title("PSD plot " + str(method))
     ax.set_xlabel('Pore width [nm]')
-    ax.set_ylabel('Distribution [dV/dw]')
+    ax.set_ylabel('Distribution, dV/dw [$cm^3 g^{-1} nm^{-1}$]')
     if labelcum:
-        ax2.set_ylabel('Cumulative Vol [$cm^3 g^{-1}$]')
+        ax2.set_ylabel('Cumulative volume [$cm^3 g^{-1}$]')
 
     lns = l1
     if labelcum:
diff --git a/tests/characterisation/test_psd_meso.py b/tests/characterisation/test_psd_meso.py
index 5ecd92bd..cbda497a 100644
--- a/tests/characterisation/test_psd_meso.py
+++ b/tests/characterisation/test_psd_meso.py
@@ -54,7 +54,7 @@ def test_psd_meso_checks(self, basic_pointisotherm):
         'BJH',
         'DH',
     ])
-    @pytest.mark.parametrize('sample', [sample for sample in DATA])
+    @pytest.mark.parametrize('sample', list(DATA))
     def test_psd_meso(self, sample, method):
         """Test psd calculation with several model isotherms."""
         sample = DATA[sample]
@@ -64,13 +64,10 @@ def test_psd_meso(self, sample, method):
             filepath = DATA_N77_PATH / sample['file']
             isotherm = pgp.isotherm_from_json(filepath)
 
-            result_dict = pmes.psd_mesoporous(
-                isotherm, psd_model=method, branch='des'
-            )
+            result_dict = pmes.psd_mesoporous(isotherm, psd_model=method, branch='des')
 
             loc = np.where(
-                result_dict['pore_distribution'] ==
-                max(result_dict['pore_distribution'])
+                result_dict['pore_distribution'] == max(result_dict['pore_distribution'])
             )
             principal_peak = result_dict['pore_widths'][loc]
 
@@ -78,8 +75,7 @@ def test_psd_meso(self, sample, method):
             err_absolute = 0.01  # 0.01
 
             assert np.isclose(
-                principal_peak, sample['psd_meso_pore_size'], err_relative,
-                err_absolute
+                principal_peak, sample['psd_meso_pore_size'], err_relative, err_absolute
             )
 
     @cleanup

From 6e76e74d2a5a18ade08bd70131befb674eaf4aa8 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Mon, 27 Jun 2022 18:16:30 +0200
Subject: [PATCH 05/10] chore: isort

---
 src/pygaps/characterisation/psd_meso.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pygaps/characterisation/psd_meso.py b/src/pygaps/characterisation/psd_meso.py
index aa1fbda1..f8facf58 100644
--- a/src/pygaps/characterisation/psd_meso.py
+++ b/src/pygaps/characterisation/psd_meso.py
@@ -6,6 +6,7 @@
 
 import numpy
 
+from pygaps import logger
 from pygaps.characterisation.models_kelvin import get_kelvin_model
 from pygaps.characterisation.models_kelvin import get_meniscus_geometry
 from pygaps.characterisation.models_thickness import get_thickness_model
@@ -15,7 +16,6 @@
 from pygaps.utilities.exceptions import CalculationError
 from pygaps.utilities.exceptions import ParameterError
 from pygaps.utilities.exceptions import pgError
-from pygaps import logger
 
 _MESO_PSD_MODELS = ['pygaps-DH', 'BJH', 'DH']
 _PORE_GEOMETRIES = ['slit', 'cylinder', 'halfopen-cylinder', 'sphere']

From 7f80aea440707acd9cde47e0de402256d1b1f4a1 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Wed, 29 Jun 2022 10:23:44 +0200
Subject: [PATCH 06/10] chore: keep flake8 happy

---
 src/pygaps/characterisation/psd_meso.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pygaps/characterisation/psd_meso.py b/src/pygaps/characterisation/psd_meso.py
index f8facf58..2ef74d3e 100644
--- a/src/pygaps/characterisation/psd_meso.py
+++ b/src/pygaps/characterisation/psd_meso.py
@@ -763,7 +763,7 @@ def psd_dollimore_heal(
         # Calculate the volume desorbed from thinning of all pores previously emptied
         # dt * \sum A - t * dt * 2 * \pi * \sum L
         d_thickness_volume = d_thickness[i] * sum_area_factor - \
-            d_thickness[i] * avg_thickness[i] * sum_2pi_length_factor # [cm3/mat]
+            d_thickness[i] * avg_thickness[i] * sum_2pi_length_factor  # [cm3/mat]
 
         # Pore volume, then store
         # dVp = (dV - dVt) * Rp

From 94b1507eb9d4c6b66ffb019d0049386c7f40d7a3 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Wed, 29 Jun 2022 22:03:24 +0200
Subject: [PATCH 07/10] fix: optional argument

---
 src/pygaps/cli/cli.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/pygaps/cli/cli.py b/src/pygaps/cli/cli.py
index 12886bc5..63336674 100644
--- a/src/pygaps/cli/cli.py
+++ b/src/pygaps/cli/cli.py
@@ -33,6 +33,7 @@ def main():
     prs.add_argument(
         'iso',
         metavar='isotherm',
+        nargs='?',
         type=pathlib.Path,
         help='isotherm to display or process',
     )

From a62204fbfd5be67bc87a4de0801bedb16e76964e Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Wed, 29 Jun 2022 22:03:32 +0200
Subject: [PATCH 08/10] chore: changelog

---
 CHANGELOG.rst | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index a5188741..204fe1e9 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -2,6 +2,12 @@
 Changelog
 =========
 
+4.2.0 (2022-06-29)
+------------------
+
+* Better error handling of nonphysical mesoporous PSD models
+* CLI now prints package version with `--version` argument
+
 4.1.1 (2022-03-10)
 ------------------
 

From 37548bf075f5b666ab65ce58809c9b405c3be911 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Wed, 29 Jun 2022 22:05:00 +0200
Subject: [PATCH 09/10] =?UTF-8?q?'4.1.2=20=E2=86=92=204.2.0'?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 conda_recipe/meta.yaml |  2 +-
 pyproject.toml         |  2 +-
 setup.cfg              | 20 ++++++++++----------
 src/pygaps/__init__.py |  2 +-
 4 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/conda_recipe/meta.yaml b/conda_recipe/meta.yaml
index bb2c6615..03169819 100644
--- a/conda_recipe/meta.yaml
+++ b/conda_recipe/meta.yaml
@@ -1,5 +1,5 @@
 {% set name = "pygaps" %}
-{% set version = "4.1.2" %}
+{% set version = "4.2.0" %}
 
 package:
   name: {{ name|lower }}
diff --git a/pyproject.toml b/pyproject.toml
index 600f74c2..030a67fd 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,7 +10,7 @@ build-backend = "setuptools.build_meta"
 [tool.setuptools_scm]
 write_to = "src/pygaps/_version.py"
 local_scheme = 'dirty-tag'
-fallback_version = '4.1.2'
+fallback_version = '4.2.0'
 
 # linting: pylint
 [tool.pylint.basic]
diff --git a/setup.cfg b/setup.cfg
index 15729f87..dbf52ab5 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 4.1.2
+current_version = 4.2.0
 commit = True
 tag = False
 message = '{current_version} → {new_version}'
@@ -16,10 +16,10 @@ author_email = mail@pauliacomi.com
 license = MIT license
 license_files = LICENSE
 url = https://github.com/pauliacomi/pygaps
-project_urls =
+project_urls = 
 	Documentation = https://pygaps.readthedocs.io
 	Source Code = https://github.com/pauliacomi/pygaps
-classifiers =
+classifiers = 
 	Development Status :: 5 - Production/Stable
 	License :: OSI Approved :: MIT License
 	Operating System :: OS Independent
@@ -38,18 +38,18 @@ classifiers =
 	Topic :: Scientific/Engineering :: Chemistry
 
 [options]
-package_dir =
+package_dir = 
 	=src
 packages = find:
 python_requires = >=3.6
 zip_safe = True
 include_package_data = True
-tests_require =
+tests_require = 
 	pytest
 	pytest-cov
 	coverage[toml]
 	nose
-install_requires =
+install_requires = 
 	numpy >= 1.16.5
 	scipy >= 1.4
 	pandas
@@ -71,11 +71,11 @@ pygaps.data = *
 pygaps.data.kernels = *.csv
 
 [options.entry_points]
-console_scripts =
+console_scripts = 
 	pygaps = pygaps.cli:main
 
 [options.extras_require]
-dev =
+dev = 
 	pytest
 	pytest-cov
 	coverage[toml]
@@ -89,7 +89,7 @@ dev =
 	bump2version
 	pygments
 	docutils
-docs =
+docs = 
 	pygments >= 2.11
 	docutils >= 0.11
 	doc8
@@ -118,7 +118,7 @@ search = version = "{current_version}"
 replace = version = "{new_version}"
 
 [flake8]
-ignore =
+ignore = 
 	E402
 	W503
 	W504
diff --git a/src/pygaps/__init__.py b/src/pygaps/__init__.py
index 4eee5fd3..3375116b 100644
--- a/src/pygaps/__init__.py
+++ b/src/pygaps/__init__.py
@@ -9,7 +9,7 @@
     from ._version import version
     __version__ = version
 except ImportError:
-    __version__ = '4.1.2'
+    __version__ = '4.2.0'
 
 import sys
 from .logging import logger

From b6f217e09347f4dcfd1eeaa2c98d453cc6039b86 Mon Sep 17 00:00:00 2001
From: Paul Iacomi <pauliacomi@users.noreply.github.com>
Date: Wed, 29 Jun 2022 22:06:23 +0200
Subject: [PATCH 10/10] ci: version bump

---
 setup.cfg | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index dbf52ab5..f76b6ef8 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -6,7 +6,7 @@ message = '{current_version} → {new_version}'
 
 [metadata]
 name = pygaps
-version = 4.1.2
+version = 4.2.0
 description = A framework for processing adsorption data for porous materials.
 long_description = file: README.rst
 long_description_content_type = text/x-rst
@@ -16,10 +16,10 @@ author_email = mail@pauliacomi.com
 license = MIT license
 license_files = LICENSE
 url = https://github.com/pauliacomi/pygaps
-project_urls = 
+project_urls =
 	Documentation = https://pygaps.readthedocs.io
 	Source Code = https://github.com/pauliacomi/pygaps
-classifiers = 
+classifiers =
 	Development Status :: 5 - Production/Stable
 	License :: OSI Approved :: MIT License
 	Operating System :: OS Independent
@@ -38,18 +38,18 @@ classifiers =
 	Topic :: Scientific/Engineering :: Chemistry
 
 [options]
-package_dir = 
+package_dir =
 	=src
 packages = find:
 python_requires = >=3.6
 zip_safe = True
 include_package_data = True
-tests_require = 
+tests_require =
 	pytest
 	pytest-cov
 	coverage[toml]
 	nose
-install_requires = 
+install_requires =
 	numpy >= 1.16.5
 	scipy >= 1.4
 	pandas
@@ -71,11 +71,11 @@ pygaps.data = *
 pygaps.data.kernels = *.csv
 
 [options.entry_points]
-console_scripts = 
+console_scripts =
 	pygaps = pygaps.cli:main
 
 [options.extras_require]
-dev = 
+dev =
 	pytest
 	pytest-cov
 	coverage[toml]
@@ -89,7 +89,7 @@ dev =
 	bump2version
 	pygments
 	docutils
-docs = 
+docs =
 	pygments >= 2.11
 	docutils >= 0.11
 	doc8
@@ -118,7 +118,7 @@ search = version = "{current_version}"
 replace = version = "{new_version}"
 
 [flake8]
-ignore = 
+ignore =
 	E402
 	W503
 	W504