From a3fbc94eff4221fd27bec8d803fa7d8e5f380a03 Mon Sep 17 00:00:00 2001
From: Martin Vonk <66305055+martinvonk@users.noreply.github.com>
Date: Fri, 10 Mar 2023 16:00:10 +0100
Subject: [PATCH 01/21] Update README.md

add badges
---
 README.md | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/README.md b/README.md
index 6a938b1..d71fef8 100644
--- a/README.md
+++ b/README.md
@@ -2,6 +2,19 @@
 
 *from Greek: πέδον, pedon -> soil*
 
+[![PyPI](https://img.shields.io/pypi/v/pedon?style=flat-square)](https://pypi.org/project/pedon/)
+[![PyPi Supported Python Versions](https://img.shields.io/pypi/pyversions/pedon?style=flat-square)](https://pypi.org/project/pedon/)
+[![Code Size](https://img.shields.io/github/languages/code-size/martinvonk/pedon?style=flat-square)](https://pypi.org/project/pedon/)
+[![PyPi Downloads](https://img.shields.io/pypi/dm/pedon?style=flat-square) ![License](https://img.shields.io/pypi/l/pedon?style=flat-square)](https://pypi.org/project/pedon/)
+
+[![Tests](https://img.shields.io/github/actions/workflow/status/martinvonk/pedon/tests.yaml?style=flat-square)](https://github.com/martinvonk/pedon/actions/workflows/tests.yaml)
+[![MyPy](https://img.shields.io/badge/type_checker-mypy-2A6DB2?style=flat-square)](https://mypy-lang.org/)
+[![Format: isort](https://img.shields.io/badge/imports-isort-ef8336?style=flat-square)](https://pycqa.github.io/isort/index.html)
+[![Format: Black](https://img.shields.io/badge/code_style-black-black?style=flat-square)](https://github.com/psf/black)
+[![Linter: flake8](https://img.shields.io/badge/linter-flake8-yellowgreen?style=flat-square)](https://flake8.pycqa.org/)
+[![Linter: ruff](https://img.shields.io/badge/linter-ruff-red?style=flat-square)](https://github.com/charliermarsh/ruff)
+
+
 Python package for (unsaturated) soil properties including pedotransfer functions. This package takes an object-oriented approach to soils, soil samples and soil models. Soil models that are available in this package are:
   - [Mualem-Van Genuchten](http://www.soilphysics.okstate.edu/teaching/soil-6583/references-folder/van%20Genuchten%201980.pdf)
   - [Brooks-Corey](https://www.wipp.energy.gov/library/cra/2009_cra/references/Others/Brooks_Corey_1964_Hydraulic_Properties_ERMS241117.pdf)

From 7eb3e6adc8e7a39ae790883d83381bc5317d44c1 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Wed, 15 Mar 2023 11:59:27 +0100
Subject: [PATCH 02/21] update Panday with sy and repr update

---
 src/pedon/soilmodel.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index 2ff1716..13fe269 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -1,4 +1,4 @@
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 from typing import Protocol
 
 import matplotlib.pyplot as plt
@@ -150,15 +150,19 @@ class Panday:
     """
 
     k_s: float
-    theta_r: float
-    theta_s: float
+    theta_r: float = field(repr=False)
+    theta_s: float = field(repr=False)
     alpha: float  # alpha
     beta: float  # n
     brook: float  # brooks-corey l
+    sr: float = field(init=False)
+    gamma: float = field(init=False, repr=False)  # 1 - 1 / beta
+    sy: float = field(init=False, repr=False)
 
     def __post_init__(self):
-        self.gamma = 1 - 1 / self.beta  # m
         self.sr = self.theta_r / self.theta_s  # theta_r / theta_s
+        self.gamma = 1 - 1 / self.beta  # m
+        self.sy = self.theta_s - self.theta_r
 
     def theta(self, h: FloatArray) -> FloatArray:
         return (self.sr + self.s(h) * (1 - self.sr)) * self.theta_s

From 6b963eac20cc13859d27db6ab7b0aa58e41ba5f7 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Wed, 15 Mar 2023 13:17:11 +0100
Subject: [PATCH 03/21] update definition of specific yield mfusg

specific yield = effective porosity - field capacity porosity
effective porosity = theta_s - theta_r
---
 src/pedon/soilmodel.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index 13fe269..48bc07f 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -162,7 +162,7 @@ class Panday:
     def __post_init__(self):
         self.sr = self.theta_r / self.theta_s  # theta_r / theta_s
         self.gamma = 1 - 1 / self.beta  # m
-        self.sy = self.theta_s - self.theta_r
+        self.sy = self.theta_s - self.theta_r - self.theta(10**2)
 
     def theta(self, h: FloatArray) -> FloatArray:
         return (self.sr + self.s(h) * (1 - self.sr)) * self.theta_s

From 6ed4493411492599820533f3119a14230804b109 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Wed, 15 Mar 2023 13:26:58 +0100
Subject: [PATCH 04/21] add function that list all possible soil names

---
 src/pedon/soil.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/pedon/soil.py b/src/pedon/soil.py
index 8e98506..69409f5 100644
--- a/src/pedon/soil.py
+++ b/src/pedon/soil.py
@@ -421,6 +421,12 @@ def from_name(self, sm: Type[SoilModel]) -> "Soil":
         self.__setattr__("model", sm(**ser))
         return self
 
+    @staticmethod
+    def list_names(sm: Type[SoilModel]) -> list[str]:
+        path = Path(__file__).parent / "datasets/Soil_Parameters.xlsx"
+        names = read_excel(path, sheet_name=sm.__name__).loc[:, "name"].to_list()
+        return names
+
     def from_staring(self, year: str = "2018") -> "Soil":
         path = Path(__file__).parent / f"datasets/Staring_{year}.xlsx"
         parameters = read_excel(path, sheet_name="parameters", index_col=0)

From 982b81ccf57ecf89e1577f9212f4d5200c4d116e Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Wed, 15 Mar 2023 13:35:47 +0100
Subject: [PATCH 05/21] make soilmodel class runtime checkable

---
 src/pedon/soilmodel.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index 48bc07f..b79c699 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -1,5 +1,5 @@
 from dataclasses import dataclass, field
-from typing import Protocol
+from typing import Protocol, runtime_checkable
 
 import matplotlib.pyplot as plt
 from numpy import abs as npabs
@@ -7,7 +7,7 @@
 
 from ._typing import FloatArray
 
-
+@runtime_checkable
 class SoilModel(Protocol):
     def theta(self, h: FloatArray) -> FloatArray:
         """Method to calculate the soil moisture content from the pressure head h"""

From c607d204d91390850d107871585813c77471adb8 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Wed, 15 Mar 2023 13:51:22 +0100
Subject: [PATCH 06/21] build in some checks

---
 src/pedon/soil.py      | 42 +++++++++++++++++++++++++++++++++---------
 src/pedon/soilmodel.py | 14 +++++++++++++-
 2 files changed, 46 insertions(+), 10 deletions(-)

diff --git a/src/pedon/soil.py b/src/pedon/soil.py
index 69409f5..1162bbd 100644
--- a/src/pedon/soil.py
+++ b/src/pedon/soil.py
@@ -9,7 +9,7 @@
 
 from ._params import get_params
 from ._typing import FloatArray
-from .soilmodel import Brooks, Genuchten, SoilModel
+from .soilmodel import Brooks, Genuchten, SoilModel, get_soilmodel
 
 
 @dataclass
@@ -406,13 +406,23 @@ class Soil:
     source: str | None = None
     description: str | None = None
 
-    def from_name(self, sm: Type[SoilModel]) -> "Soil":
+    def from_name(self, sm: Type[SoilModel] | SoilModel | str) -> "Soil":
+        if isinstance(sm, SoilModel):
+            smn = sm.__class__.name
+            sm = type(sm)
+        elif isinstance(sm, str):
+            smn = sm
+            sm = get_soilmodel(smn)
+        elif issubclass(sm, SoilModel):
+            smn = sm.__name__
+        else:
+            raise ValueError(
+                f"Argument must either be Type[SoilModel] | SoilModel | str,"
+                f"not {type(sm)}"
+            )
+
         path = Path(__file__).parent / "datasets/Soil_Parameters.xlsx"
-        ser = (
-            read_excel(path, sheet_name=sm.__name__, index_col=0)
-            .loc[self.name]
-            .to_dict()
-        )
+        ser = read_excel(path, sheet_name=smn, index_col=0).loc[self.name].to_dict()
         # path = Path(__file__).parent / f"datasets/{sm.__name__}.csv"
         # ser = read_csv(path, index_col=["name"]).loc[self.name].to_dict()
         self.__setattr__("type", ser.pop("soil type"))
@@ -422,12 +432,26 @@ def from_name(self, sm: Type[SoilModel]) -> "Soil":
         return self
 
     @staticmethod
-    def list_names(sm: Type[SoilModel]) -> list[str]:
+    def list_names(sm: Type[SoilModel] | SoilModel | str) -> list[str]:
+        if isinstance(sm, SoilModel):
+            smn = sm.__class__.name
+        elif isinstance(sm, str):
+            smn = sm
+        elif issubclass(sm, SoilModel):
+            smn = sm.__name__
+        else:
+            raise ValueError(
+                f"Argument must either be Type[SoilModel] | SoilModel | str,"
+                f"not {type(sm)}"
+            )
+
         path = Path(__file__).parent / "datasets/Soil_Parameters.xlsx"
-        names = read_excel(path, sheet_name=sm.__name__).loc[:, "name"].to_list()
+        names = read_excel(path, sheet_name=smn).loc[:, "name"].to_list()
         return names
 
     def from_staring(self, year: str = "2018") -> "Soil":
+        if year not in ("2001", "2018"):
+            raise ValueError(f"Year must either be '2001' or '2018, not {year}")
         path = Path(__file__).parent / f"datasets/Staring_{year}.xlsx"
         parameters = read_excel(path, sheet_name="parameters", index_col=0)
         ser = parameters.loc[self.name].to_dict()
diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index b79c699..70c2449 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -1,5 +1,5 @@
 from dataclasses import dataclass, field
-from typing import Protocol, runtime_checkable
+from typing import Protocol, runtime_checkable, Type
 
 import matplotlib.pyplot as plt
 from numpy import abs as npabs
@@ -7,6 +7,7 @@
 
 from ._typing import FloatArray
 
+
 @runtime_checkable
 class SoilModel(Protocol):
     def theta(self, h: FloatArray) -> FloatArray:
@@ -243,6 +244,17 @@ def plot(self):
         return plot_swrc(self)
 
 
+def get_soilmodel(soilmodel_name: str) -> Type[SoilModel]:
+    sms = {
+        "Genuchten": Genuchten,
+        "Brooks": Brooks,
+        "Gardner": Gardner,
+        "Panday": Panday,
+        "Fredlund": Fredlund,
+    }
+    return sms[soilmodel_name]
+
+
 def plot_swrc(
     sm: SoilModel,
     saturation: bool = False,

From e59a5bb3d873f201dcef444dfdc436a75cdb55a7 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 16 Mar 2023 11:32:05 +0100
Subject: [PATCH 07/21] fix isinstance test

runtime_checkable does not differentiate between a class and instance for the ifinstance check. so now check on attribute __name__ which is only an attribute for a class
---
 src/pedon/soil.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/src/pedon/soil.py b/src/pedon/soil.py
index 1162bbd..cf61dcf 100644
--- a/src/pedon/soil.py
+++ b/src/pedon/soil.py
@@ -408,13 +408,15 @@ class Soil:
 
     def from_name(self, sm: Type[SoilModel] | SoilModel | str) -> "Soil":
         if isinstance(sm, SoilModel):
-            smn = sm.__class__.name
-            sm = type(sm)
+            if hasattr(sm, "__name__"):
+                smn = sm.__name__
+            else:
+                smn = sm.__class__.__name__
+                sm = type(sm)
         elif isinstance(sm, str):
             smn = sm
             sm = get_soilmodel(smn)
-        elif issubclass(sm, SoilModel):
-            smn = sm.__name__
+
         else:
             raise ValueError(
                 f"Argument must either be Type[SoilModel] | SoilModel | str,"
@@ -434,11 +436,15 @@ def from_name(self, sm: Type[SoilModel] | SoilModel | str) -> "Soil":
     @staticmethod
     def list_names(sm: Type[SoilModel] | SoilModel | str) -> list[str]:
         if isinstance(sm, SoilModel):
-            smn = sm.__class__.name
+            if hasattr(sm, "__name__"):
+                smn = sm.__name__
+            else:
+                smn = sm.__class__.__name__
+                sm = type(sm)
         elif isinstance(sm, str):
             smn = sm
-        elif issubclass(sm, SoilModel):
-            smn = sm.__name__
+            sm = get_soilmodel(smn)
+
         else:
             raise ValueError(
                 f"Argument must either be Type[SoilModel] | SoilModel | str,"

From d5f6aeca44db22ec36f0951ad50c37375c5e8b45 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 16 Mar 2023 15:25:42 +0100
Subject: [PATCH 08/21] parse ax to SoilModel.plot method

---
 src/pedon/soilmodel.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index 70c2449..de6631c 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -22,7 +22,7 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         """Method to calcualte the permeability from the pressure head h"""
         ...
 
-    def plot(self):
+    def plot(self, ax: plt.Axes | None = None):
         """Method to plot the soil water retention curve"""
         ...
 
@@ -61,8 +61,8 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             s = self.s(h)
         return self.k_s * s**self.l * (1 - (1 - s ** (1 / self.m)) ** self.m) ** 2
 
-    def plot(self):
-        return plot_swrc(self)
+    def plot(self, ax: plt.Axes | None = None):
+        return plot_swrc(self, ax=ax)
 
 
 @dataclass
@@ -109,8 +109,8 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             s = self.s(h)
         return self.k_s * s ** (3 + 2 / self.l)
 
-    def plot(self):
-        return plot_swrc(self)
+    def plot(self, ax: plt.Axes | None = None):
+        return plot_swrc(self, ax=ax)
 
 
 @dataclass
@@ -139,8 +139,8 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             return self.k_s * self.a * theta**self.m
         return self.k_s * (self.a / (self.b + npabs(h) ** self.m))
 
-    def plot(self):
-        return plot_swrc(self)
+    def plot(self, ax: plt.Axes | None = None):
+        return plot_swrc(self, ax=ax)
 
 
 @dataclass
@@ -176,8 +176,8 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             s = self.s(h)
         return self.k_s * s**self.brook
 
-    def plot(self):
-        return plot_swrc(self)
+    def plot(self, ax: plt.Axes | None = None):
+        return plot_swrc(self, ax=ax)
 
 
 @dataclass
@@ -240,8 +240,8 @@ def theta_d(
 
         return self.k_s * (teller / noemer)
 
-    def plot(self):
-        return plot_swrc(self)
+    def plot(self, ax: plt.Axes | None = None):
+        return plot_swrc(self, ax=ax)
 
 
 def get_soilmodel(soilmodel_name: str) -> Type[SoilModel]:

From edc0bcec59af04f022e954697d7a6d61f1a6f710 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 23 Mar 2023 11:39:49 +0100
Subject: [PATCH 09/21] add bubblept to panday

---
 src/pedon/soilmodel.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index de6631c..c1766df 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -156,7 +156,8 @@ class Panday:
     alpha: float  # alpha
     beta: float  # n
     brook: float  # brooks-corey l
-    sr: float = field(init=False)
+    h_b: float = field(default=0, repr=False)
+    sr: float = field(init=False, repr=True)
     gamma: float = field(init=False, repr=False)  # 1 - 1 / beta
     sy: float = field(init=False, repr=False)
 
@@ -169,7 +170,7 @@ def theta(self, h: FloatArray) -> FloatArray:
         return (self.sr + self.s(h) * (1 - self.sr)) * self.theta_s
 
     def s(self, h: FloatArray) -> FloatArray:
-        return (1 + npabs(self.alpha * h) ** self.beta) ** -self.gamma
+        return (1 + npabs(self.alpha * (h - self.h_b)) ** self.beta) ** -self.gamma
 
     def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         if s is None:

From cb9dc9000b239ce1af0059d18d76a1d3456d2d05 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 23 Mar 2023 16:22:56 +0100
Subject: [PATCH 10/21] fix protocol since plot returns axes, not None

---
 src/pedon/soilmodel.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index c1766df..8c438ad 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -1,5 +1,5 @@
 from dataclasses import dataclass, field
-from typing import Protocol, runtime_checkable, Type
+from typing import Protocol, Type, runtime_checkable
 
 import matplotlib.pyplot as plt
 from numpy import abs as npabs
@@ -22,7 +22,7 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         """Method to calcualte the permeability from the pressure head h"""
         ...
 
-    def plot(self, ax: plt.Axes | None = None):
+    def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         """Method to plot the soil water retention curve"""
         ...
 
@@ -61,7 +61,7 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             s = self.s(h)
         return self.k_s * s**self.l * (1 - (1 - s ** (1 / self.m)) ** self.m) ** 2
 
-    def plot(self, ax: plt.Axes | None = None):
+    def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
 
 
@@ -109,7 +109,7 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             s = self.s(h)
         return self.k_s * s ** (3 + 2 / self.l)
 
-    def plot(self, ax: plt.Axes | None = None):
+    def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
 
 
@@ -139,7 +139,7 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             return self.k_s * self.a * theta**self.m
         return self.k_s * (self.a / (self.b + npabs(h) ** self.m))
 
-    def plot(self, ax: plt.Axes | None = None):
+    def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
 
 
@@ -177,7 +177,7 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
             s = self.s(h)
         return self.k_s * s**self.brook
 
-    def plot(self, ax: plt.Axes | None = None):
+    def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
 
 
@@ -241,7 +241,7 @@ def theta_d(
 
         return self.k_s * (teller / noemer)
 
-    def plot(self, ax: plt.Axes | None = None):
+    def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
 
 

From 5138492bcf7f96481e361915f54b3b66d9740455 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 23 Mar 2023 16:24:07 +0100
Subject: [PATCH 11/21] Update __init__.py

https://github.com/microsoft/pylance-release/issues/2953#issuecomment-1168408943
---
 src/pedon/__init__.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/pedon/__init__.py b/src/pedon/__init__.py
index 3677d6f..5bc4911 100644
--- a/src/pedon/__init__.py
+++ b/src/pedon/__init__.py
@@ -2,5 +2,9 @@
 # flake8: noqa
 from . import _params, soil, soilmodel
 from ._version import __version__
-from .soil import Soil, SoilSample
-from .soilmodel import Brooks, Fredlund, Gardner, Genuchten, Panday, plot_hcf, plot_swrc
+from .soil import Soil as Soil
+from .soil import SoilSample as SoilSample
+from .soilmodel import Brooks as Brooks
+from .soilmodel import Genuchten as Genuchten
+from .soilmodel import Panday as Panday
+from .soilmodel import plot_hcf, plot_swrc

From d8ac9d1a45c83474052eb2d9c105291b71c22df3 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Fri, 24 Mar 2023 16:15:51 +0100
Subject: [PATCH 12/21] improve brook parameter initial settings

---
 src/pedon/_params.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pedon/_params.py b/src/pedon/_params.py
index e61abaa..cb11fd1 100644
--- a/src/pedon/_params.py
+++ b/src/pedon/_params.py
@@ -74,7 +74,7 @@
             "theta_s": 0.4,
             "alpha": 0.02,
             "beta": 2.3,
-            "brook": 0.5,
+            "brook": 10.0,
         },
         "p_min": {
             "k_s": 0.001,
@@ -82,7 +82,7 @@
             "theta_s": 0.2,
             "alpha": 0.001,
             "beta": 1.0,
-            "brook": -7,
+            "brook": 1.0,
         },
         "p_max": {
             "k_s": 100000.0,
@@ -90,7 +90,7 @@
             "theta_s": 0.5,
             "alpha": 0.15,
             "beta": 12,
-            "brook": 8,
+            "brook": 40.0,
         },
         "swrc": {
             "k_s": False,

From a80bbf40394390d2e2b7c969cefe6b7ff7c7cef5 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Fri, 24 Mar 2023 19:46:32 +0100
Subject: [PATCH 13/21] Update _params.py

better max params for panday
---
 src/pedon/_params.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pedon/_params.py b/src/pedon/_params.py
index cb11fd1..cd1eefb 100644
--- a/src/pedon/_params.py
+++ b/src/pedon/_params.py
@@ -22,7 +22,7 @@
             "k_s": 100000.0,
             "theta_r": 0.2,
             "theta_s": 0.9,
-            "alpha": 0.15,
+            "alpha": 0.20,
             "n": 12,
             "l": 8,
         },
@@ -88,9 +88,9 @@
             "k_s": 100000.0,
             "theta_r": 0.2,
             "theta_s": 0.5,
-            "alpha": 0.15,
+            "alpha": 0.30,
             "beta": 12,
-            "brook": 40.0,
+            "brook": 50.0,
         },
         "swrc": {
             "k_s": False,

From 54da345e698737490a647a599873de3b90daa57d Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 30 Mar 2023 16:23:53 +0200
Subject: [PATCH 14/21] add relative permeability functions

---
 src/pedon/soilmodel.py | 40 +++++++++++++++++++++++++++++-----------
 1 file changed, 29 insertions(+), 11 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index 8c438ad..439dba9 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -18,6 +18,10 @@ def s(self, h: FloatArray) -> FloatArray:
         """Method to calculate the effective saturation from the pressure head h"""
         ...
 
+    def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+        """Method to calcualte the relative permeability from the pressure head h"""
+        ...
+
     def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         """Method to calcualte the permeability from the pressure head h"""
         ...
@@ -56,10 +60,13 @@ def theta(self, h: FloatArray) -> FloatArray:
     def s(self, h: FloatArray) -> FloatArray:
         return (self.theta(h) - self.theta_r) / (self.theta_s - self.theta_r)
 
-    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+    def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         if s is None:
             s = self.s(h)
-        return self.k_s * s**self.l * (1 - (1 - s ** (1 / self.m)) ** self.m) ** 2
+        return s**self.l * (1 - (1 - s ** (1 / self.m)) ** self.m) ** 2
+
+    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+        return self.k_s * self.k_r(h=h, s=s)
 
     def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
@@ -104,10 +111,13 @@ def s(self, h: FloatArray) -> FloatArray:
             s[h >= self.h_b] = (h[h >= self.h_b] / self.h_b) ** -self.l
             return s
 
-    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+    def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         if s is None:
             s = self.s(h)
-        return self.k_s * s ** (3 + 2 / self.l)
+        return s ** (3 + 2 / self.l)
+
+    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+        return self.k_s * self.k_r(h=h, s=s)
 
     def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
@@ -133,11 +143,14 @@ def theta(self, h: FloatArray) -> FloatArray:
     def s(self, h: FloatArray) -> FloatArray:
         return (self.theta(h) - self.theta_r) / (self.theta_s - self.theta_r)
 
-    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+    def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         if s is not None:
             theta = s * (self.theta_s - self.theta_r) + self.theta_r
-            return self.k_s * self.a * theta**self.m
-        return self.k_s * (self.a / (self.b + npabs(h) ** self.m))
+            return self.a * theta**self.m
+        return self.a / (self.b + npabs(h) ** self.m)
+
+    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+        return self.k_s * self.k_r(h=h, s=s)
 
     def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
@@ -172,10 +185,13 @@ def theta(self, h: FloatArray) -> FloatArray:
     def s(self, h: FloatArray) -> FloatArray:
         return (1 + npabs(self.alpha * (h - self.h_b)) ** self.beta) ** -self.gamma
 
-    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+    def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
         if s is None:
             s = self.s(h)
-        return self.k_s * s**self.brook
+        return s**self.brook
+
+    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+        return self.k_s * self.k_r(h=h, s=s)
 
     def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)
@@ -201,7 +217,7 @@ def theta(self, h: FloatArray) -> FloatArray:
     def s(self, h: FloatArray) -> FloatArray:
         return self.theta(h) / self.theta_s
 
-    def k(self, h: FloatArray, s: FloatArray | None = None):
+    def k_r(self, h: FloatArray, s: FloatArray | None = None):
         if s is not None:
             raise NotImplementedError(
                 "Can only calculate the hydraulic conductivity"
@@ -238,8 +254,10 @@ def theta_d(
                 / exp(y)
                 * theta_d(exp(y), self.a, self.n, self.m, self.theta_s)
             )
+        return teller / noemer
 
-        return self.k_s * (teller / noemer)
+    def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray:
+        return self.k_s * self.k_r(h=h, s=s)
 
     def plot(self, ax: plt.Axes | None = None) -> plt.Axes:
         return plot_swrc(self, ax=ax)

From 2a2b3191d1e28efd006a42f6d90cf258ed125e7a Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 30 Mar 2023 18:33:11 +0200
Subject: [PATCH 15/21] plot relative_k method

---
 src/pedon/soilmodel.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index 439dba9..f868d72 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -303,6 +303,7 @@ def plot_swrc(
 def plot_hcf(
     sm: SoilModel,
     ax: plt.Axes | None = None,
+    relative_k: bool = False,
     **kwargs: dict,
 ) -> plt.Axes:
     """Plot the hydraulic conductivity function"""
@@ -314,7 +315,10 @@ def plot_hcf(
 
     h = logspace(-6, 10, num=1000)
 
-    k = sm.k(h=h)
+    if relative_k:
+        k = sm.k_r(h=h)
+    else:
+        k = sm.k(h=h)
 
     ax.plot(k, h, label=sm.__class__.__name__, **kwargs)
     ax.set_ylim(1e-3, 1e6)

From 48bc40ae55a4125e1222e09fa3d8e032ded2db92 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 30 Mar 2023 18:35:13 +0200
Subject: [PATCH 16/21] add method to fit on relative permeability

---
 src/pedon/soil.py | 40 ++++++++++++++++++++++++++++++----------
 1 file changed, 30 insertions(+), 10 deletions(-)

diff --git a/src/pedon/soil.py b/src/pedon/soil.py
index cf61dcf..5bcaf05 100644
--- a/src/pedon/soil.py
+++ b/src/pedon/soil.py
@@ -116,14 +116,26 @@ def fit(
         W1: float | None = None,
         W2: float | None = None,
         return_res: bool = False,
+        relative_k: bool = False,
     ) -> SoilModel:
         """Same method as RETC"""
+
+        theta = self.theta
+        k = self.k
+
         if pbounds is None:
             pbounds = get_params(sm.__name__)
-            pbounds.loc["k_s", "p_ini"] = max(self.k)
-            pbounds.loc["k_s", "p_max"] = max(self.k) * 10
-            pbounds.loc["theta_s", "p_ini"] = max(self.theta)
-            pbounds.loc["theta_s", "p_max"] = max(self.theta) + 0.02
+            if relative_k:
+                pbounds = pbounds.drop("k_s")
+            else:
+                pbounds.loc["k_s", "p_ini"] = max(k)
+                pbounds.loc["k_s", "p_max"] = max(k) * 10
+            pbounds.loc["theta_s", "p_ini"] = max(theta)
+            pbounds.loc["theta_s", "p_max"] = max(theta) + 0.02
+
+        if relative_k:
+            k_s = max(k)
+            k /= k_s
 
         if weights is None:
             weights = ones(self.h.shape)
@@ -132,14 +144,20 @@ def fit(
             W1 = 0.1
 
         if W2 is None:
-            M = len(self.k) + len(self.theta)
-            N = len(self.theta)
-            W2 = (M - N) * sum(weights * self.theta) / (N * sum(weights * self.k))
+            M = len(k) + len(theta)
+            N = len(theta)
+            W2 = (M - N) * sum(weights * theta) / (N * sum(weights * k))
 
         def fit_staring(p: FloatArray) -> FloatArray:
-            sml = sm(**dict(zip(pbounds.index, p)))
-            theta_diff = sml.theta(h=self.h) - self.theta
-            k_diff = log(sml.k(h=self.h)) - log(self.k)
+            est_pars = dict(zip(pbounds.index, p))
+            if relative_k:
+                est_pars["k_s"] = k_s
+            sml = sm(**est_pars)
+            theta_diff = sml.theta(h=self.h) - theta
+            if relative_k:
+                k_diff = log(sml.k_r(h=self.h)) - log(k)
+            else:
+                k_diff = log(sml.k(h=self.h)) - log(k)
             diff = append(weights * theta_diff, weights * W1 * W2 * k_diff)
             return diff
 
@@ -152,6 +170,8 @@ def fit_staring(p: FloatArray) -> FloatArray:
             ),
         )
         opt_pars = dict(zip(pbounds.index, res.x))
+        if relative_k:
+            opt_pars["k_s"] = k_s
         opt_sm = sm(**opt_pars)
         if return_res:
             return opt_sm, {"res": res}

From 879f8abbca0764f2d3659c517e1df55d5a9cd572 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 30 Mar 2023 18:51:49 +0200
Subject: [PATCH 17/21] remove relative_k plotting option

---
 src/pedon/soilmodel.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index f868d72..a7135c6 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -303,7 +303,6 @@ def plot_swrc(
 def plot_hcf(
     sm: SoilModel,
     ax: plt.Axes | None = None,
-    relative_k: bool = False,
     **kwargs: dict,
 ) -> plt.Axes:
     """Plot the hydraulic conductivity function"""
@@ -314,11 +313,7 @@ def plot_hcf(
         ax.set_xscale("log")
 
     h = logspace(-6, 10, num=1000)
-
-    if relative_k:
-        k = sm.k_r(h=h)
-    else:
-        k = sm.k(h=h)
+    k = sm.k(h=h)
 
     ax.plot(k, h, label=sm.__class__.__name__, **kwargs)
     ax.set_ylim(1e-3, 1e6)

From 5a4ce8257f786ce84a4f6f97429b8d3d07863d5b Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Thu, 30 Mar 2023 19:09:47 +0200
Subject: [PATCH 18/21] make optional method for providing k_s

---
 src/pedon/soil.py | 17 +++++------------
 1 file changed, 5 insertions(+), 12 deletions(-)

diff --git a/src/pedon/soil.py b/src/pedon/soil.py
index 5bcaf05..194c584 100644
--- a/src/pedon/soil.py
+++ b/src/pedon/soil.py
@@ -116,7 +116,7 @@ def fit(
         W1: float | None = None,
         W2: float | None = None,
         return_res: bool = False,
-        relative_k: bool = False,
+        k_s: float | None = None,
     ) -> SoilModel:
         """Same method as RETC"""
 
@@ -125,7 +125,7 @@ def fit(
 
         if pbounds is None:
             pbounds = get_params(sm.__name__)
-            if relative_k:
+            if k_s is not None:
                 pbounds = pbounds.drop("k_s")
             else:
                 pbounds.loc["k_s", "p_ini"] = max(k)
@@ -133,10 +133,6 @@ def fit(
             pbounds.loc["theta_s", "p_ini"] = max(theta)
             pbounds.loc["theta_s", "p_max"] = max(theta) + 0.02
 
-        if relative_k:
-            k_s = max(k)
-            k /= k_s
-
         if weights is None:
             weights = ones(self.h.shape)
 
@@ -150,14 +146,11 @@ def fit(
 
         def fit_staring(p: FloatArray) -> FloatArray:
             est_pars = dict(zip(pbounds.index, p))
-            if relative_k:
+            if k_s is not None:
                 est_pars["k_s"] = k_s
             sml = sm(**est_pars)
             theta_diff = sml.theta(h=self.h) - theta
-            if relative_k:
-                k_diff = log(sml.k_r(h=self.h)) - log(k)
-            else:
-                k_diff = log(sml.k(h=self.h)) - log(k)
+            k_diff = log(sml.k(h=self.h)) - log(k)
             diff = append(weights * theta_diff, weights * W1 * W2 * k_diff)
             return diff
 
@@ -170,7 +163,7 @@ def fit_staring(p: FloatArray) -> FloatArray:
             ),
         )
         opt_pars = dict(zip(pbounds.index, res.x))
-        if relative_k:
+        if k_s is not None:
             opt_pars["k_s"] = k_s
         opt_sm = sm(**opt_pars)
         if return_res:

From e96a985ec73503edb72b482021d07329986fb4f2 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Tue, 11 Apr 2023 14:39:40 +0200
Subject: [PATCH 19/21] Update soilmodel.py

make sure float is actually float
---
 src/pedon/soilmodel.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py
index a7135c6..dc4ec18 100644
--- a/src/pedon/soilmodel.py
+++ b/src/pedon/soilmodel.py
@@ -169,7 +169,7 @@ class Panday:
     alpha: float  # alpha
     beta: float  # n
     brook: float  # brooks-corey l
-    h_b: float = field(default=0, repr=False)
+    h_b: float = field(default=0.0, repr=False)
     sr: float = field(init=False, repr=True)
     gamma: float = field(init=False, repr=False)  # 1 - 1 / beta
     sy: float = field(init=False, repr=False)

From b55e6513c007e4141647aa0ff1eab7f41aa8d676 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Mon, 15 May 2023 08:46:06 +0200
Subject: [PATCH 20/21] bump version

---
 src/pedon/_version.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pedon/_version.py b/src/pedon/_version.py
index 3b93d0b..27fdca4 100644
--- a/src/pedon/_version.py
+++ b/src/pedon/_version.py
@@ -1 +1 @@
-__version__ = "0.0.2"
+__version__ = "0.0.3"

From c9da775326c1b31504be67f6e17074d3228ca032 Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Mon, 15 May 2023 08:46:17 +0200
Subject: [PATCH 21/21] better year check

---
 src/pedon/soil.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pedon/soil.py b/src/pedon/soil.py
index 194c584..b2828ed 100644
--- a/src/pedon/soil.py
+++ b/src/pedon/soil.py
@@ -469,8 +469,8 @@ def list_names(sm: Type[SoilModel] | SoilModel | str) -> list[str]:
         return names
 
     def from_staring(self, year: str = "2018") -> "Soil":
-        if year not in ("2001", "2018"):
-            raise ValueError(f"Year must either be '2001' or '2018, not {year}")
+        if year not in ("2001", 2001, "2018", 2018):
+            raise ValueError(f"Year must either be '2001' or '2018', not {year}")
         path = Path(__file__).parent / f"datasets/Staring_{year}.xlsx"
         parameters = read_excel(path, sheet_name="parameters", index_col=0)
         ser = parameters.loc[self.name].to_dict()