photosynthesis.py

Go to the documentation of this file.

"""
Photosynthesis parameter structures and data classes for PyHelios.
 
This module provides Python data structures that mirror the C++ parameter
classes used by the PhotosynthesisModel plugin, with proper defaults and
validation support.
"""
 
from dataclasses import dataclass, field
from typing import List, Optional, Union
import math
 
# Known species in the photosynthesis library (21 species with aliases)
PHOTOSYNTHESIS_SPECIES = [
    "Almond", "Apple", "Cherry", "Prune", "Pear", 
    "PistachioFemale", "PistachioMale", "Walnut",
    "Grape",  # cv. Cabernet Sauvignon
    "Elderberry", "Toyon", "Big_Leaf_Maple", "Western_Redbud", "Baylaurel", "Olive",
    "EasternRedbudSunlit", "EasternRedbudShaded"
]
 
# Species aliases for case-insensitive and format-flexible lookup
SPECIES_ALIASES = {
    # Standard names (already in PHOTOSYNTHESIS_SPECIES)
    "almond": "Almond",
    "apple": "Apple", 
    "cherry": "Cherry",
    "prune": "Prune",
    "pear": "Pear",
    "pistachiofemale": "PistachioFemale",
    "pistachiomale": "PistachioMale",
    "walnut": "Walnut",
    "grape": "Grape",
    "elderberry": "Elderberry",
    "toyon": "Toyon",
    "big_leaf_maple": "Big_Leaf_Maple",
    "western_redbud": "Western_Redbud",
    "baylaurel": "Baylaurel",
    "olive": "Olive",
    "easternredbudsunlit": "EasternRedbudSunlit",
    "easternredbudshaded": "EasternRedbudShaded",
    
    # Common aliases and variations
    "bigleafmaple": "Big_Leaf_Maple",
    "bigmaple": "Big_Leaf_Maple",
    "westernredbud": "Western_Redbud",
    "redbud": "Western_Redbud",
    "easternredbud": "EasternRedbudSunlit",  # Default to sunlit
    "pistachio": "PistachioFemale",  # Default to female
    "cabernet": "Grape",
    "cabernetSauvignon": "Grape",
    "grapevine": "Grape"
}
 
 
# Topt sentinel: PhotosyntheticTemperatureResponseParameters stores Topt in Kelvin and uses
# 10000 K to mean "no optimum". After conversion to Celsius that is ~9726.85, far above the
# 100 C that helios-core's own validateOptimalTemperature accepts, so any value at or above
# 200 C in the flat array must be the sentinel rather than user data. Mirrors
# C4_NO_OPTIMUM_TOPT_C in native/src/pyhelios_wrapper_photosynthesis.cpp.
_NO_OPTIMUM_TOPT_C = 200.0
 
 
@dataclass
class PhotosyntheticTemperatureResponseParameters:
    """
    Temperature response parameters for photosynthetic processes.
    
    These parameters define how photosynthetic rates vary with temperature
    using the modified Arrhenius equation.
    
    Attributes:
        value_at_25C: Value of the parameter at 25°C
        dHa: Activation energy (rate of increase parameter)
        dHd: Deactivation energy (rate of decrease parameter)  
        Topt: Optimum temperature in Kelvin (10000K means no optimum)
    """
    value_at_25C: float = 100.0
    dHa: float = 60.0
    dHd: float = 600.0
    Topt: float = 10000.0  # Very high = no temperature optimum
    
    def __post_init__(self):
        """Validate parameter values after initialization."""
        if not math.isfinite(self.value_at_25C):
            raise ValueError("value_at_25C must be finite")
        if not math.isfinite(self.dHa) or self.dHa < 0:
            raise ValueError("dHa must be finite and non-negative")
        if not math.isfinite(self.dHd) or self.dHd < 0:
            raise ValueError("dHd must be finite and non-negative")
        if not math.isfinite(self.Topt) or self.Topt < 0:
            raise ValueError("Topt must be finite and non-negative")
        # Mirrors validateDeactivationEnergy in helios-core 1.3.80: the peaked Arrhenius form
        # evaluates ln(dHd/dHa - 1), which is undefined unless dHd > dHa. dHa <= 0 means no
        # Arrhenius term is applied, so dHd is unused.
        if self.dHa > 0 and self.dHd <= self.dHa:
            raise ValueError(
                f"Deactivation energy dHd must be strictly greater than activation energy dHa "
                f"for a peaked temperature response. Received dHa = {self.dHa} kJ/mol and "
                f"dHd = {self.dHd} kJ/mol. The peaked Arrhenius form evaluates "
                f"ln(dHd/dHa - 1), which is undefined when dHd <= dHa. Increase dHd "
                f"(a typical value is 10*dHa, or 200-600 kJ/mol) or omit dHd to use the default."
            )
 
 
@dataclass
class EmpiricalModelCoefficients:
    """
    Empirical photosynthesis model coefficients.
    
    This model uses empirical relationships to estimate photosynthetic
    rates based on environmental conditions.
    
    Attributes:
        Tref: Reference temperature (K)
        Ci_ref: Reference CO2 concentration (μmol CO2/mol air)
        Asat: Light-saturated photosynthetic rate (μmol/m²/s)
        theta: Half-saturation light level (W/m²)
        Tmin: Minimum temperature for photosynthesis (K)
        Topt: Optimum temperature for photosynthesis (K)
        q: Temperature response parameter (unitless)
        R: Respiration temperature coefficient (μmol·K^0.5/m²/s)
        ER: Respiration activation energy (1/K)
        kC: CO2 response coefficient (unitless)
    """
    Tref: float = 298.0      # K
    Ci_ref: float = 290.0    # μmol CO2/mol air
    Asat: float = 18.18      # μmol/m²/s
    theta: float = 62.03     # W/m²
    Tmin: float = 290.0      # K
    Topt: float = 303.0      # K  
    q: float = 0.344         # unitless
    R: float = 1.663e5       # μmol·K^0.5/m²/s
    ER: float = 3740.0       # 1/K
    kC: float = 0.791        # unitless
    
    def __post_init__(self):
        """Validate parameter values after initialization."""
        if self.Tref <= 0:
            raise ValueError("Reference temperature must be positive")
        if self.Ci_ref <= 0:
            raise ValueError("Reference CO2 concentration must be positive")
        if self.Asat < 0:
            raise ValueError("Light-saturated photosynthetic rate cannot be negative")
        if self.theta <= 0:
            raise ValueError("Half-saturation light level must be positive")
        if self.Tmin <= 0:
            raise ValueError("Minimum temperature must be positive")
        if self.Topt <= 0:
            raise ValueError("Optimum temperature must be positive")
        if self.Tmin >= self.Topt:
            raise ValueError("Minimum temperature must be less than optimum temperature")
        if self.q <= 0:
            raise ValueError("Temperature response parameter must be positive")
        if self.R < 0:
            raise ValueError("Respiration coefficient cannot be negative")
        if self.ER < 0:
            raise ValueError("Respiration activation energy cannot be negative")
        if self.kC < 0:
            raise ValueError("CO2 response coefficient cannot be negative")
        # Mirrors the coefficient checks helios-core 1.3.80 added to the empirical
        # temperature response f_T. Before 1.3.80 the Tmin/Topt/Tref/q coefficients were
        # completely inert; now they drive f_T, and these two combinations make its
        # reference denominator zero, giving an infinite or NaN assimilation rate.
        if self.Tref <= self.Tmin:
            raise ValueError(
                f"Reference temperature Tref ({self.Tref} K) must be greater than the minimum "
                f"temperature Tmin ({self.Tmin} K); otherwise the empirical temperature "
                f"response f_T has a zero or negative reference denominator."
            )
        denom_ref = (1.0 + self.q) * self.Topt - self.Tmin - self.q * self.Tref
        if abs(denom_ref) < 1e-6:
            raise ValueError(
                f"Empirical temperature response coefficients are degenerate: "
                f"(1+q)*Topt - Tmin - q*Tref = {denom_ref}, which makes the reference "
                f"denominator of f_T zero. Adjust Topt, Tmin, Tref or q."
            )
    
    def to_array(self) -> List[float]:
        """Convert to float array for C++ interface."""
        return [
            self.Tref, self.Ci_ref, self.Asat, self.theta, self.Tmin,
            self.Topt, self.q, self.R, self.ER, self.kC
        ]
    
    @classmethod
    def from_array(cls, coefficients: List[float]) -> 'EmpiricalModelCoefficients':
        """Create from float array (from C++ interface)."""
        if len(coefficients) < 10:
            raise ValueError("Need at least 10 coefficients for empirical model")
        return cls(
            Tref=coefficients[0], Ci_ref=coefficients[1], Asat=coefficients[2],
            theta=coefficients[3], Tmin=coefficients[4], Topt=coefficients[5],
            q=coefficients[6], R=coefficients[7], ER=coefficients[8], kC=coefficients[9]
        )
 
 
@dataclass  
class FarquharModelCoefficients:
    """
    Farquhar-von Caemmerer-Berry photosynthesis model coefficients.
    
    This model provides a mechanistic description of leaf photosynthesis
    based on biochemical limitations and temperature responses.
    
    Core Parameters (at 25°C):
        Vcmax: Maximum carboxylation rate (μmol/m²/s, -1 = uninitialized)
        Jmax: Maximum electron transport rate (μmol/m²/s, -1 = uninitialized)  
        alpha: Quantum efficiency of photosystem II (μmol electrons/μmol photons)
        Rd: Dark respiration rate (μmol/m²/s, -1 = uninitialized)
        O: Ambient oxygen concentration (mmol/mol)
        TPU_flag: Enable triose phosphate utilization limitation (0/1)
        
    Temperature Response Parameters:
        c_*: Scaling factor for Arrhenius equation
        dH_*: Activation energy for temperature response
    """
    # Core parameters at 25°C
    Vcmax: float = -1.0      # μmol/m²/s (uninitialized)
    Jmax: float = -1.0       # μmol/m²/s (uninitialized)
    alpha: float = -1.0      # unitless (uninitialized)
    Rd: float = -1.0         # μmol/m²/s (uninitialized)
    O: float = 213.5         # ambient oxygen concentration (mmol/mol)
    TPU_flag: int = 0        # enable TPU limitation
    
    # Temperature scaling factors (c_*)
    c_Rd: float = 18.72
    c_Vcmax: float = 26.35
    c_Jmax: float = 18.86
    c_Gamma: float = 19.02
    c_Kc: float = 38.05
    c_Ko: float = 20.30
    
    # Activation energies (dH_*)  
    dH_Rd: float = 46.39
    dH_Vcmax: float = 65.33
    dH_Jmax: float = 46.36
    dH_Gamma: float = 37.83
    dH_Kc: float = 79.43
    dH_Ko: float = 36.38
 
    # Mesophyll conductance gm (helios-core 1.3.72+). The default math.inf reproduces
    # the legacy Cc = Ci behaviour (no mesophyll diffusion limitation). When packed into
    # the flat coefficient array, slots [18..21] are (gm_at_25C, dHa, Topt_C, dHd) using
    # the same -1 sentinel convention as the C4 model: dHa < 0 → constant gm with no
    # temperature response.
    gm_at_25C: float = float('inf')   # mol CO2 / m^2 / s / bar
    dHa_gm: float = -1.0              # kJ/mol; -1 disables temperature response
    Topt_gm_C: float = -1.0           # Celsius; -1 → monotonic Arrhenius
    dHd_gm: float = -1.0              # kJ/mol; -1 → default
 
    # Temperature response parameter containers
    _vcmax_temp_response: Optional[PhotosyntheticTemperatureResponseParameters] = field(default=None, init=False)
    _jmax_temp_response: Optional[PhotosyntheticTemperatureResponseParameters] = field(default=None, init=False)
    _rd_temp_response: Optional[PhotosyntheticTemperatureResponseParameters] = field(default=None, init=False)
    _alpha_temp_response: Optional[PhotosyntheticTemperatureResponseParameters] = field(default=None, init=False)
    _theta_temp_response: Optional[PhotosyntheticTemperatureResponseParameters] = field(default=None, init=False)
    
    def __post_init__(self):
        """Validate parameter values after initialization."""
        if self.O <= 0:
            raise ValueError("Oxygen concentration must be positive")
        if self.TPU_flag not in (0, 1):
            raise ValueError("TPU_flag must be 0 or 1")
        
        # Validate temperature parameters
        for param_name, value in [
            ('c_Rd', self.c_Rd), ('c_Vcmax', self.c_Vcmax), ('c_Jmax', self.c_Jmax),
            ('c_Gamma', self.c_Gamma), ('c_Kc', self.c_Kc), ('c_Ko', self.c_Ko),
            ('dH_Rd', self.dH_Rd), ('dH_Vcmax', self.dH_Vcmax), ('dH_Jmax', self.dH_Jmax),
            ('dH_Gamma', self.dH_Gamma), ('dH_Kc', self.dH_Kc), ('dH_Ko', self.dH_Ko)
        ]:
            if not math.isfinite(value):
                raise ValueError(f"Temperature parameter {param_name} must be finite")
    
    def setVcmax(self, vcmax_at_25c: float, dha: Optional[float] = None, 
                 topt: Optional[float] = None, dhd: Optional[float] = None) -> None:
        """Set Vcmax with temperature response (mimics C++ overloads)."""
        if dha is None:
            # 1-parameter version
            self._vcmax_temp_response = PhotosyntheticTemperatureResponseParameters(vcmax_at_25c)
        elif topt is None:
            # 2-parameter version
            self._vcmax_temp_response = PhotosyntheticTemperatureResponseParameters(vcmax_at_25c, dha)
        elif dhd is None:
            # 3-parameter version
            # C++ defaults dHd to 10*dHa for the 3-argument peaked form.
            self._vcmax_temp_response = PhotosyntheticTemperatureResponseParameters(
                vcmax_at_25c, dHa=dha, dHd=10.0 * dha, Topt=273.15 + topt)
        else:
            # 4-parameter version
            self._vcmax_temp_response = PhotosyntheticTemperatureResponseParameters(
                vcmax_at_25c, dHa=dha, dHd=dhd, Topt=273.15 + topt)
        
        self.Vcmax = vcmax_at_25c
    
    def setJmax(self, jmax_at_25c: float, dha: Optional[float] = None,
                topt: Optional[float] = None, dhd: Optional[float] = None) -> None:
        """Set Jmax with temperature response (mimics C++ overloads)."""
        if dha is None:
            self._jmax_temp_response = PhotosyntheticTemperatureResponseParameters(jmax_at_25c)
        elif topt is None:
            self._jmax_temp_response = PhotosyntheticTemperatureResponseParameters(jmax_at_25c, dha)
        elif dhd is None:
            # C++ defaults dHd to 10*dHa for the 3-argument peaked form.
            self._jmax_temp_response = PhotosyntheticTemperatureResponseParameters(
                jmax_at_25c, dHa=dha, dHd=10.0 * dha, Topt=273.15 + topt)
        else:
            self._jmax_temp_response = PhotosyntheticTemperatureResponseParameters(
                jmax_at_25c, dHa=dha, dHd=dhd, Topt=273.15 + topt)
        
        self.Jmax = jmax_at_25c
    
    def setRd(self, rd_at_25c: float, dha: Optional[float] = None,
              topt: Optional[float] = None, dhd: Optional[float] = None) -> None:
        """Set dark respiration with temperature response (mimics C++ overloads)."""
        if dha is None:
            self._rd_temp_response = PhotosyntheticTemperatureResponseParameters(rd_at_25c)
        elif topt is None:
            self._rd_temp_response = PhotosyntheticTemperatureResponseParameters(rd_at_25c, dha)
        elif dhd is None:
            # C++ defaults dHd to 10*dHa for the 3-argument peaked form.
            self._rd_temp_response = PhotosyntheticTemperatureResponseParameters(
                rd_at_25c, dHa=dha, dHd=10.0 * dha, Topt=273.15 + topt)
        else:
            self._rd_temp_response = PhotosyntheticTemperatureResponseParameters(
                rd_at_25c, dHa=dha, dHd=dhd, Topt=273.15 + topt)
        
        self.Rd = rd_at_25c
    
    def setQuantumEfficiency_alpha(self, alpha_at_25c: float, dha: Optional[float] = None,
                                   topt: Optional[float] = None, dhd: Optional[float] = None) -> None:
        """Set quantum efficiency with temperature response (mimics C++ overloads).""" 
        if dha is None:
            self._alpha_temp_response = PhotosyntheticTemperatureResponseParameters(alpha_at_25c)
        elif topt is None:
            self._alpha_temp_response = PhotosyntheticTemperatureResponseParameters(alpha_at_25c, dha)
        elif dhd is None:
            # C++ defaults dHd to 10*dHa for the 3-argument peaked form.
            self._alpha_temp_response = PhotosyntheticTemperatureResponseParameters(
                alpha_at_25c, dHa=dha, dHd=10.0 * dha, Topt=273.15 + topt)
        else:
            self._alpha_temp_response = PhotosyntheticTemperatureResponseParameters(
                alpha_at_25c, dHa=dha, dHd=dhd, Topt=273.15 + topt)
        
        self.alpha = alpha_at_25c
    
    def setLightResponseCurvature_theta(self, theta_at_25c: float, dha: Optional[float] = None,
                                        topt: Optional[float] = None, dhd: Optional[float] = None) -> None:
        """Set light response curvature with temperature response (mimics C++ overloads)."""
        if dha is None:
            self._theta_temp_response = PhotosyntheticTemperatureResponseParameters(theta_at_25c)
        elif topt is None:
            self._theta_temp_response = PhotosyntheticTemperatureResponseParameters(theta_at_25c, dha)
        elif dhd is None:
            # C++ defaults dHd to 10*dHa for the 3-argument peaked form.
            self._theta_temp_response = PhotosyntheticTemperatureResponseParameters(
                theta_at_25c, dHa=dha, dHd=10.0 * dha, Topt=273.15 + topt)
        else:
            self._theta_temp_response = PhotosyntheticTemperatureResponseParameters(
                theta_at_25c, dHa=dha, dHd=dhd, Topt=273.15 + topt)
    
    def getVcmaxTempResponse(self) -> PhotosyntheticTemperatureResponseParameters:
        """Get Vcmax temperature response parameters."""
        if self._vcmax_temp_response is None:
            return PhotosyntheticTemperatureResponseParameters(self.Vcmax if self.Vcmax > 0 else 100.0)
        return self._vcmax_temp_response
    
    def getJmaxTempResponse(self) -> PhotosyntheticTemperatureResponseParameters:
        """Get Jmax temperature response parameters."""
        if self._jmax_temp_response is None:
            return PhotosyntheticTemperatureResponseParameters(self.Jmax if self.Jmax > 0 else 100.0)
        return self._jmax_temp_response
    
    def getRdTempResponse(self) -> PhotosyntheticTemperatureResponseParameters:
        """Get dark respiration temperature response parameters."""
        if self._rd_temp_response is None:
            return PhotosyntheticTemperatureResponseParameters(self.Rd if self.Rd > 0 else 1.0)
        return self._rd_temp_response
    
    def getQuantumEfficiencyTempResponse(self) -> PhotosyntheticTemperatureResponseParameters:
        """Get quantum efficiency temperature response parameters."""
        if self._alpha_temp_response is None:
            return PhotosyntheticTemperatureResponseParameters(self.alpha if self.alpha > 0 else 0.3)
        return self._alpha_temp_response
    
    def getLightResponseCurvatureTempResponse(self) -> PhotosyntheticTemperatureResponseParameters:
        """Get light response curvature temperature response parameters.""" 
        if self._theta_temp_response is None:
            return PhotosyntheticTemperatureResponseParameters(0.7)
        return self._theta_temp_response
    
    def _temp_response_block(self, response: PhotosyntheticTemperatureResponseParameters,
                             fallback_value: float) -> List[float]:
        """Pack a temperature response into its 4-float (value, dHa, Topt_C, dHd) block.
 
        ``Topt`` is stored in Kelvin and defaults to 10000 K to mean "no optimum". The flat
        array carries Topt in Celsius with -1 as the no-optimum sentinel, matching
        ``packTempResponse`` in native/src/pyhelios_wrapper_photosynthesis.cpp.
        """
        if response is None:
            return [fallback_value, -1.0, -1.0, -1.0]
        topt_c = response.Topt - 273.15
        if topt_c >= _NO_OPTIMUM_TOPT_C:
            topt_c = -1.0
        return [response.value_at_25C, response.dHa, topt_c, response.dHd]
 
    def to_array(self) -> List[float]:
        """Convert to the flat float array used by the native interface.
 
        Slots 0..17 are the legacy Farquhar fields (Vcmax/Jmax/alpha/Rd/O/TPU_flag plus
        the 12 c_*/dH_* temperature constants). Slots 18..21 carry the mesophyll
        conductance gm temperature response: (gm_at_25C, dHa, Topt_C, dHd) using the
        -1 sentinel convention (dHa < 0 → constant, Topt_C < 0 → monotonic Arrhenius,
        dHd < 0 → default deactivation energy). Slots 22..37 carry the same 4-float block
        for Vcmax, Jmax, Rd and alpha in that order, and slots 38..41 the same block for the
        light response curvature theta.
        """
        return [
            self.Vcmax, self.Jmax, self.alpha, self.Rd, self.O, float(self.TPU_flag),
            # Temperature scaling factors
            self.c_Vcmax, self.dH_Vcmax, self.c_Jmax, self.dH_Jmax,
            self.c_Rd, self.dH_Rd, self.c_Kc, self.dH_Kc,
            self.c_Ko, self.dH_Ko, self.c_Gamma, self.dH_Gamma,
            # Mesophyll conductance gm temperature response (1.3.72+)
            self.gm_at_25C, self.dHa_gm, self.Topt_gm_C, self.dHd_gm,
            # Full rate temperature responses (1.3.80+)
            *self._temp_response_block(self._vcmax_temp_response, self.Vcmax),
            *self._temp_response_block(self._jmax_temp_response, self.Jmax),
            *self._temp_response_block(self._rd_temp_response, self.Rd),
            *self._temp_response_block(self._alpha_temp_response, self.alpha),
            # Light response curvature theta temperature response (slots 38..41). Theta has no
            # legacy scalar slot, so this block is the only representation that crosses.
            *self._temp_response_block(self._theta_temp_response, 0.7),
        ]
 
    @classmethod
    def from_array(cls, coefficients: List[float]) -> 'FarquharModelCoefficients':
        """Create from float array (from C++ interface).
 
        Accepts the legacy 18-float layout (pre-1.3.72), the 22-float layout with mesophyll
        conductance gm in slots 18..21, and the 38-float layout (1.3.80+) that additionally
        carries the Vcmax/Jmax/Rd/alpha temperature responses in slots 22..37, and the
        42-float layout that adds the light response curvature theta in slots 38..41. Shorter
        arrays leave the corresponding responses unset, reproducing the earlier behaviour.
        """
        if len(coefficients) < 18:
            raise ValueError("Need at least 18 coefficients for Farquhar model")
 
        gm_at_25C = coefficients[18] if len(coefficients) > 18 else float('inf')
        dHa_gm = coefficients[19] if len(coefficients) > 19 else -1.0
        Topt_gm_C = coefficients[20] if len(coefficients) > 20 else -1.0
        dHd_gm = coefficients[21] if len(coefficients) > 21 else -1.0
 
        instance = cls(
            Vcmax=coefficients[0], Jmax=coefficients[1], alpha=coefficients[2],
            Rd=coefficients[3], O=coefficients[4], TPU_flag=int(coefficients[5]),
            c_Vcmax=coefficients[6], dH_Vcmax=coefficients[7],
            c_Jmax=coefficients[8], dH_Jmax=coefficients[9],
            c_Rd=coefficients[10], dH_Rd=coefficients[11],
            c_Kc=coefficients[12], dH_Kc=coefficients[13],
            c_Ko=coefficients[14], dH_Ko=coefficients[15],
            c_Gamma=coefficients[16], dH_Gamma=coefficients[17],
            gm_at_25C=gm_at_25C, dHa_gm=dHa_gm, Topt_gm_C=Topt_gm_C, dHd_gm=dHd_gm,
        )
 
        if len(coefficients) >= 38:
            for offset, setter in ((22, instance.setVcmax), (26, instance.setJmax),
                                   (30, instance.setRd), (34, instance.setQuantumEfficiency_alpha)):
                value, dha, topt, dhd = coefficients[offset:offset + 4]
                if dha < 0:
                    setter(value)
                elif topt < 0:
                    setter(value, dha)
                elif dhd < 0:
                    setter(value, dha, topt)
                else:
                    setter(value, dha, topt, dhd)
 
        if len(coefficients) >= 42:
            value, dha, topt, dhd = coefficients[38:42]
            if dha < 0:
                instance.setLightResponseCurvature_theta(value)
            elif topt < 0:
                instance.setLightResponseCurvature_theta(value, dha)
            elif dhd < 0:
                instance.setLightResponseCurvature_theta(value, dha, topt)
            else:
                instance.setLightResponseCurvature_theta(value, dha, topt, dhd)
 
        return instance
 
 
def validate_species_name(species: str) -> str:
    """
    Validate and normalize species name for photosynthesis library.
    
    Args:
        species: Species name (case insensitive, supports aliases)
        
    Returns:
        Normalized species name
        
    Raises:
        ValueError: If species is not recognized
    """
    if not species:
        raise ValueError("Species name cannot be empty")
    
    # Try exact match first (case sensitive)
    if species in PHOTOSYNTHESIS_SPECIES:
        return species
    
    # Try case-insensitive match
    species_lower = species.lower()
    if species_lower in SPECIES_ALIASES:
        return SPECIES_ALIASES[species_lower]
    
    # Try case-insensitive match against known species
    for known_species in PHOTOSYNTHESIS_SPECIES:
        if known_species.lower() == species_lower:
            return known_species
    
    # Species not found - provide helpful error message
    available_species = sorted(set(list(PHOTOSYNTHESIS_SPECIES) + list(SPECIES_ALIASES.keys())))
    raise ValueError(
        f"Unknown species '{species}'. Available species and aliases:\n"
        f"  {', '.join(available_species[:8])}\n"
        f"  {', '.join(available_species[8:16])}\n"
        f"  {', '.join(available_species[16:])}"
    )
 
 
def get_available_species() -> List[str]:
    """Get list of available species in the photosynthesis library."""
    return sorted(PHOTOSYNTHESIS_SPECIES.copy())
 
 
def get_species_aliases() -> dict:
    """Get dictionary of species aliases."""
    return SPECIES_ALIASES.copy()