PhotosynthesisModel.h

Go to the documentation of this file.

 
#ifndef PHOTOSYNTHESIS_MODEL
#define PHOTOSYNTHESIS_MODEL
 
#include "Context.h"
 
#include <limits>
 
struct PhotosyntheticTemperatureResponseParameters {
private:
    static void validateOptimalTemperature(float optimum_temperature_in_C) {
        if (optimum_temperature_in_C < 0.f) {
            helios::helios_runtime_error("ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = " + std::to_string(optimum_temperature_in_C) +
                                         " C. Please check that temperature is provided in units of Celsius, not Kelvin.");
        }
        if (optimum_temperature_in_C > 100.f) {
            helios::helios_runtime_error("ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = " + std::to_string(optimum_temperature_in_C) +
                                         " C. This value is biologically unrealistic and likely indicates temperature was provided in Kelvin instead of Celsius. Please convert to Celsius (subtract 273.15 from Kelvin value).");
        }
    }
 
 
    static void validateDeactivationEnergy(float rate_of_increase_dHa, float rate_of_decrease_dHd) {
        if (rate_of_increase_dHa <= 0.f) {
            return; // No Arrhenius term is applied, so dHd is unused.
        }
        if (rate_of_decrease_dHd <= rate_of_increase_dHa) {
            helios::helios_runtime_error("ERROR (PhotosyntheticTemperatureResponseParameters): Deactivation energy dHd must be strictly greater than activation energy dHa for a peaked temperature response. Received dHa = " +
                                         std::to_string(rate_of_increase_dHa) + " kJ/mol and dHd = " + std::to_string(rate_of_decrease_dHd) +
                                         " kJ/mol. The peaked Arrhenius form evaluates ln(dHd/dHa - 1), which is undefined when dHd <= dHa. Increase dHd (a typical value is 10*dHa, or 200-600 kJ/mol) or omit dHd to use the default.");
        }
    }
 
public:
    PhotosyntheticTemperatureResponseParameters() {
        value_at_25C = 100.0f;
        dHa = 60.0f;
        dHd = 600.0f;
        Topt = 10000.f;
    }
 
    PhotosyntheticTemperatureResponseParameters(float value_at_25C) {
        this->value_at_25C = value_at_25C;
        this->dHa = 0.f;
        this->dHd = 600.f;
        this->Topt = 10000.f;
    }
 
    PhotosyntheticTemperatureResponseParameters(float value_at_25C, float rate_of_increase_dHa) {
        this->value_at_25C = value_at_25C;
        this->dHa = rate_of_increase_dHa;
        if (rate_of_increase_dHa > 0.f) {
            this->dHd = 10.f * rate_of_increase_dHa;
        } else {
            this->dHd = 600.f;
        }
        this->Topt = 10000.f;
    }
 
    PhotosyntheticTemperatureResponseParameters(float value_at_25C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        validateOptimalTemperature(optimum_temperature_in_C);
        this->value_at_25C = value_at_25C;
        this->dHa = rate_of_increase_dHa;
        if (rate_of_increase_dHa > 0.f) {
            this->dHd = 10.f * rate_of_increase_dHa;
        } else {
            this->dHd = 600.f;
        }
        this->Topt = 273.15f + optimum_temperature_in_C;
    }
 
    PhotosyntheticTemperatureResponseParameters(float value_at_25C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        validateOptimalTemperature(optimum_temperature_in_C);
        validateDeactivationEnergy(rate_of_increase_dHa, rate_of_decrease_dHd);
        this->value_at_25C = value_at_25C;
        this->dHa = rate_of_increase_dHa;
        this->dHd = rate_of_decrease_dHd;
        this->Topt = 273.15f + optimum_temperature_in_C;
    }
 
    float value_at_25C;
    float dHa;
    float dHd;
    float Topt;
};
 
 
struct EmpiricalModelCoefficients {
 
    EmpiricalModelCoefficients() {
        Tref = 298; // K
        Ci_ref = 290; // umol CO2/mol air
        Asat = 18.18; // umol/m^2-s
        theta = 62.03; // umol/m^2-s (PPFD; i_PAR is converted from W/m^2 before use)
        Tmin = 290; // K
        Topt = 303; // K
        q = 0.344; // unitless
        R = 1.663e5; // umol-K^0.5/m^2-s
        ER = 3740; // 1/K
        kC = 0.791; // unitless
    }
 
    // reference values
    float Tref;
    float Ci_ref;
 
    // light response coefficients
    float Asat;
    float theta;
 
    // assimilation temperature response coefficients
    float Tmin;
    float Topt;
    float q;
 
    // respiration temperature response coefficients
    float R;
    float ER;
 
    // CO2 response coefficients
    float kC;
};
 
struct FarquharModelCoefficients {
public:
    FarquharModelCoefficients() {
 
        // parameters (at TL = 25C)
        Vcmax = -1.f; // umol/m^2/s
        Jmax = -1.f; // umol/m^2/s
        alpha = -1.f; // unitless
        Rd = -1.f; // umol/m^2/s
 
        O = 213.5; // ambient oxygen concentration (mmol/mol)
 
        // temperature parameters
        c_Rd = 18.72;
        c_Vcmax = 26.35;
        c_Jmax = 18.86;
        c_Gamma = 19.02;
        c_Kc = 38.05;
        c_Ko = 20.30;
 
        dH_Rd = 46.39f;
        dH_Vcmax = 65.33f;
        dH_Jmax = 46.36f;
        dH_Gamma = 37.83;
        dH_Kc = 79.43;
        dH_Ko = 36.38;
 
        TPU_flag = 0;
    }
 
    // options
    int TPU_flag; // run model with TPU limitation
 
    // Deprecated legacy parameters. These are retained for backwards compatibility with code that
    // assigns them directly; a value > 0 selects a simple (non-peaked) Arrhenius response using the
    // corresponding c_*/dH_* coefficients below. Calling the matching setter (e.g. setVcmax())
    // resets the field to the -1 sentinel so that the setter is always authoritative — otherwise the
    // two representations could silently disagree about which value the model actually uses.
    // Prefer the setter API in new code.
    float Vcmax;
    float Jmax;
    float Rd;
    float alpha;
    float O;
 
    // temperature parameters
    float c_Rd;
    float c_Vcmax;
    float c_Jmax;
    float c_Gamma;
    float c_Kc;
    float c_Ko;
 
    float dH_Rd;
    float dH_Vcmax;
    float dH_Jmax;
    float dH_Gamma;
    float dH_Kc;
    float dH_Ko;
 
 
    void setVcmax(float Vcmax) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax);
        this->Vcmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setVcmax(float Vcmax_at_25_C, float rate_of_increase_dHa) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax_at_25_C, rate_of_increase_dHa);
        this->Vcmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setVcmax(float Vcmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        float rate_of_decrease_dHd = 10.f * rate_of_increase_dHa;
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->Vcmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setVcmax(float Vcmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->Vcmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setJmax(float Jmax) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax);
        this->Jmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setJmax(float Jmax_at_25_C, float rate_of_increase_dHa) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax_at_25_C, rate_of_increase_dHa);
        this->Jmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setJmax(float Jmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        float rate_of_decrease_dHd = 10.f * rate_of_increase_dHa;
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->Jmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setJmax(float Jmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->Jmax = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setTPU(float TPU) {
        this->TPUTempResponse = PhotosyntheticTemperatureResponseParameters(TPU);
        this->TPU_flag = 1;
    }
 
 
    void setTPU(float TPU_at_25_C, float rate_of_increase_dHa) {
        this->TPUTempResponse = PhotosyntheticTemperatureResponseParameters(TPU_at_25_C, rate_of_increase_dHa);
        this->TPU_flag = 1;
    }
 
 
    void setTPU(float TPU_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        float rate_of_decrease_dHd = 10.f * rate_of_increase_dHa;
        this->TPUTempResponse = PhotosyntheticTemperatureResponseParameters(TPU_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->TPU_flag = 1;
    }
 
 
    void setTPU(float TPU_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->TPUTempResponse = PhotosyntheticTemperatureResponseParameters(TPU_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->TPU_flag = 1;
    }
 
 
    void setRd(float Rd) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd);
        this->Rd = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setRd(float Rd_at_25_C, float rate_of_increase_dHa) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd_at_25_C, rate_of_increase_dHa);
        this->Rd = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setRd(float Rd_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        float rate_of_decrease_dHd = 10.f * rate_of_increase_dHa;
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->Rd = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setRd(float Rd_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->Rd = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setQuantumEfficiency_alpha(float alpha) {
        this->alphaTempResponse = PhotosyntheticTemperatureResponseParameters(alpha);
        this->alpha = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setQuantumEfficiency_alpha(float alpha_at_25_C, float rate_of_increase_dHa) {
        this->alphaTempResponse = PhotosyntheticTemperatureResponseParameters(alpha_at_25_C, rate_of_increase_dHa);
        this->alpha = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setQuantumEfficiency_alpha(float alpha_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        float rate_of_decrease_dHd = 10.f * rate_of_increase_dHa;
        this->alphaTempResponse = PhotosyntheticTemperatureResponseParameters(alpha_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->alpha = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setQuantumEfficiency_alpha(float alpha_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->alphaTempResponse = PhotosyntheticTemperatureResponseParameters(alpha_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
        this->alpha = -1.f; // legacy field: setter is authoritative
    }
 
 
    void setLightResponseCurvature_theta(float theta) {
        this->thetaTempResponse = PhotosyntheticTemperatureResponseParameters(theta);
    }
 
 
    void setLightResponseCurvature_theta(float theta_at_25_C, float rate_of_increase_dHa) {
        this->thetaTempResponse = PhotosyntheticTemperatureResponseParameters(theta_at_25_C, rate_of_increase_dHa);
    }
 
 
    void setLightResponseCurvature_theta(float theta_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        float rate_of_decrease_dHd = 10.f * rate_of_increase_dHa;
        this->thetaTempResponse = PhotosyntheticTemperatureResponseParameters(theta_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
 
    void setLightResponseCurvature_theta(float theta_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->thetaTempResponse = PhotosyntheticTemperatureResponseParameters(theta_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
 
    void setMesophyllConductance_gm(float gm) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm);
    }
 
 
    void setMesophyllConductance_gm(float gm_at_25_C, float rate_of_increase_dHa) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm_at_25_C, rate_of_increase_dHa);
    }
 
 
    void setMesophyllConductance_gm(float gm_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C);
    }
 
 
    void setMesophyllConductance_gm(float gm_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
    PhotosyntheticTemperatureResponseParameters getVcmaxTempResponse() {
        return this->VcmaxTempResponse;
    }
 
    PhotosyntheticTemperatureResponseParameters getJmaxTempResponse() {
        return this->JmaxTempResponse;
    }
 
    PhotosyntheticTemperatureResponseParameters getTPUTempResponse() {
        return this->TPUTempResponse;
    }
 
    PhotosyntheticTemperatureResponseParameters getRdTempResponse() {
        return this->RdTempResponse;
    }
 
    PhotosyntheticTemperatureResponseParameters getQuantumEfficiencyTempResponse() {
        return this->alphaTempResponse;
    }
 
    PhotosyntheticTemperatureResponseParameters getLightResponseCurvatureTempResponse() {
        return this->thetaTempResponse;
    }
 
    PhotosyntheticTemperatureResponseParameters getMesophyllConductance_gmTempResponse() const {
        return this->gmTempResponse;
    }
 
 
private:
    PhotosyntheticTemperatureResponseParameters VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(100);
    PhotosyntheticTemperatureResponseParameters JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(200);
    PhotosyntheticTemperatureResponseParameters TPUTempResponse = PhotosyntheticTemperatureResponseParameters(2.0);
    PhotosyntheticTemperatureResponseParameters RdTempResponse = PhotosyntheticTemperatureResponseParameters(1.5);
    PhotosyntheticTemperatureResponseParameters alphaTempResponse = PhotosyntheticTemperatureResponseParameters(0.5);
    PhotosyntheticTemperatureResponseParameters thetaTempResponse = PhotosyntheticTemperatureResponseParameters(0.f);
    PhotosyntheticTemperatureResponseParameters gmTempResponse = PhotosyntheticTemperatureResponseParameters(std::numeric_limits<float>::infinity());
 
 
    friend class PhotosynthesisModel;
};
 
 
 
struct C4ModelCoefficients {
public:
    C4ModelCoefficients() {
        // Flexible temperature responses — Setaria viridis defaults from Table 1 (Boyd et al. 2015 for Vpmax/Vcmax/Rd,
        // Ubierna et al. 2017 refit Arrhenius for gm). Jmax is re-fit from the paper's Gaussian form
        // Jmax(T) = 400·exp(−((T−43)/26)²) to the peaked-Arrhenius form used here — Jmax(25°C) = 247.69 μmol/m²/s
        // matches the spreadsheet at the 25°C reference; dHa / Topt / dHd chosen to approximate the Gaussian.
        VpmaxTempResponse = PhotosyntheticTemperatureResponseParameters(200.f, 50.1f);
        VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(40.f, 78.0f);
        JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(247.69f, 77.9f, 43.0f, 260.f);
        RdTempResponse = PhotosyntheticTemperatureResponseParameters(1.f, 66.4f);
        gmTempResponse = PhotosyntheticTemperatureResponseParameters(1.f, 49.8f);
    }
 
    // === (1) Flexible temperature-responsive parameters: setter overloads mirror C3 API ===
 
    void setVpmax(float Vpmax) {
        this->VpmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vpmax);
    }
    void setVpmax(float Vpmax_at_25_C, float rate_of_increase_dHa) {
        this->VpmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vpmax_at_25_C, rate_of_increase_dHa);
    }
    void setVpmax(float Vpmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        this->VpmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vpmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C);
    }
    void setVpmax(float Vpmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->VpmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vpmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
    void setVcmax(float Vcmax) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax);
    }
    void setVcmax(float Vcmax_at_25_C, float rate_of_increase_dHa) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax_at_25_C, rate_of_increase_dHa);
    }
    void setVcmax(float Vcmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C);
    }
    void setVcmax(float Vcmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->VcmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Vcmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
    void setJmax(float Jmax) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax);
    }
    void setJmax(float Jmax_at_25_C, float rate_of_increase_dHa) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax_at_25_C, rate_of_increase_dHa);
    }
    void setJmax(float Jmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C);
    }
    void setJmax(float Jmax_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->JmaxTempResponse = PhotosyntheticTemperatureResponseParameters(Jmax_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
    void setRd(float Rd) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd);
    }
    void setRd(float Rd_at_25_C, float rate_of_increase_dHa) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd_at_25_C, rate_of_increase_dHa);
    }
    void setRd(float Rd_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C);
    }
    void setRd(float Rd_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->RdTempResponse = PhotosyntheticTemperatureResponseParameters(Rd_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
    void setMesophyllConductance_gm(float gm) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm);
    }
    void setMesophyllConductance_gm(float gm_at_25_C, float rate_of_increase_dHa) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm_at_25_C, rate_of_increase_dHa);
    }
    void setMesophyllConductance_gm(float gm_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C);
    }
    void setMesophyllConductance_gm(float gm_at_25_C, float rate_of_increase_dHa, float optimum_temperature_in_C, float rate_of_decrease_dHd) {
        this->gmTempResponse = PhotosyntheticTemperatureResponseParameters(gm_at_25_C, rate_of_increase_dHa, optimum_temperature_in_C, rate_of_decrease_dHd);
    }
 
    // Accessors (read-only) used for material-data serialization
    PhotosyntheticTemperatureResponseParameters getVpmaxTempResponse() const { return VpmaxTempResponse; }
    PhotosyntheticTemperatureResponseParameters getVcmaxTempResponse() const { return VcmaxTempResponse; }
    PhotosyntheticTemperatureResponseParameters getJmaxTempResponse() const { return JmaxTempResponse; }
    PhotosyntheticTemperatureResponseParameters getRdTempResponse() const { return RdTempResponse; }
    PhotosyntheticTemperatureResponseParameters getMesophyllConductance_gmTempResponse() const { return gmTempResponse; }
 
    // === (2) Rubisco + PEPC kinetic constants — hardcoded 25 °C values with simple Arrhenius responses.
    //     Edit directly if different species are needed (paper defaults are for Setaria viridis). Units: μbar
    //     for Cm-like partial pressures, kJ/mol for activation energies.
 
    float Kc_25 = 1210.f;             
    float Ko_25 = 292000.f;           
    float Kp_25 = 82.f;               
    float gamma_star_25 = 3.81679e-4f;
    float Om_25 = 210000.f;           
 
    float dH_Kc = 64.2f;              
    float dH_Ko = 10.5f;              
    float dH_Kp = 38.3f;              
    // IMPORTANT — DO NOT change dH_gamma_star back to a negative value.
    // The von Caemmerer (2021) Setaria spreadsheet (C4__model_setaria__11-06-2021.xlsm, cell C17)
    // lists -31.1 kJ/mol, but this is a transcription error: Boyd, Gandin & Cousins (2015)
    // Plant Physiol 169:1850 Table I reports Ea = -31.1 kJ/mol for S_c/o (Rubisco specificity,
    // which decreases with temperature). Because γ* = 0.5/S_c/o, the activation energy for γ*
    // has the OPPOSITE sign: Ea(γ*) = -Ea(S_c/o) = +31.1 kJ/mol. The spreadsheet copied Boyd's
    // S_c/o value into the γ* row without the required sign flip. Woodford, Ermakova, Furbank
    // & von Caemmerer (2025, bioRxiv 10.1101/2025.06.03.657559) Table 1 silently corrects this
    // to +31.1 kJ/mol. Positive is also consistent with Bernacchi et al. (2001) for C3
    // (Ea(Γ*) ≈ +37.8 kJ/mol) and with the physical fact that photorespiration increases with T.
    float dH_gamma_star = 31.1f;      
    float dH_Om = 0.f;                
 
    // === (3) User-tunable scalar parameters (paper defaults in comments) ===
 
    float alpha_psII_fraction = 0.f;  
    float x_etr_partition = 0.4f;     
    float Vpr = 80.f;                 
    float Rm_frac = 0.5f;             
    float fcyc = 0.45f;               
    float gbs = 0.003f;               
    float ao = 0.047f;                
    float absorptance = 0.85f;        
    float f_spectral = 0.15f;         
    float theta_etr = 0.7f;           
    float h_protons = 4.f;            
    // Proton stoichiometries for the z factor in the Aj quadratic (ATP / linear-e⁻ ratio).
    // vC2021 originally collapsed H_J and H_Jcyc to a single "3" in the numerator of the z formula,
    // which implicitly assumes H_Jcyc = 2 (PGR5-only cyclic flow). Woodford et al. (2025) argue NDH
    // dominates Setaria cyclic flow (Ermakova 2024 methyl-viologen work), giving H_Jcyc = 3.4, and
    // derive fcyc = 0.45 consistent with Bjorkman & Demmig 1987 / Yin & Struik 2012 quantum yields.
    // Together these recover the paper Table 1 value z = 1.45 via the corrected formula in
    // computeC4RatesFromCm: z = (H_J·(1−fcyc) + H_Jcyc·fcyc) / (h_protons·(1−fcyc)).
    float H_J = 3.f;                  
    float H_Jcyc = 3.4f;              
 
private:
    PhotosyntheticTemperatureResponseParameters VpmaxTempResponse;
    PhotosyntheticTemperatureResponseParameters VcmaxTempResponse;
    PhotosyntheticTemperatureResponseParameters JmaxTempResponse;
    PhotosyntheticTemperatureResponseParameters RdTempResponse;
    PhotosyntheticTemperatureResponseParameters gmTempResponse;
 
    friend class PhotosynthesisModel;
};
 
 
class PhotosynthesisModel {
public:
 
    explicit PhotosynthesisModel(helios::Context *a_context);
 
    static int selfTest(int argc = 0, char **argv = nullptr);
 
    void setModelType_Empirical();
 
    void setModelType_Farquhar();
 
    void setModelType_C4();
 
 
    void setModelCoefficients(const EmpiricalModelCoefficients &modelcoefficients);
 
 
    void setModelCoefficients(const EmpiricalModelCoefficients &modelcoefficients, const std::vector<uint> &UUIDs);
 
 
    void setModelCoefficients(const FarquharModelCoefficients &modelcoefficients);
 
 
    void setModelCoefficients(const FarquharModelCoefficients &modelcoefficients, const std::vector<uint> &UUIDs);
 
 
    void setModelCoefficients(const std::vector<FarquharModelCoefficients> &modelcoefficients, const std::vector<uint> &UUIDs);
 
 
    void setFarquharCoefficientsFromLibrary(const std::string &species);
 
 
    void setFarquharCoefficientsFromLibrary(const std::string &species, const std::vector<uint> &UUIDs);
 
 
    FarquharModelCoefficients getFarquharCoefficientsFromLibrary(const std::string &species);
 
 
    void setModelCoefficients(const std::string &material_label, const EmpiricalModelCoefficients &coeffs);
 
 
    void setModelCoefficients(const std::string &material_label, const FarquharModelCoefficients &coeffs);
 
 
    void setFarquharCoefficientsFromLibrary(const std::string &species, const std::string &material_label);
 
 
    void setModelCoefficients(const C4ModelCoefficients &modelcoefficients);
 
 
    void setModelCoefficients(const C4ModelCoefficients &modelcoefficients, const std::vector<uint> &UUIDs);
 
 
    void setModelCoefficients(const std::string &material_label, const C4ModelCoefficients &coeffs);
 
 
    void setC4CoefficientsFromLibrary(const std::string &species);
 
 
    void setC4CoefficientsFromLibrary(const std::string &species, const std::vector<uint> &UUIDs);
 
 
    void setC4CoefficientsFromLibrary(const std::string &species, const std::string &material_label);
 
 
    C4ModelCoefficients getC4CoefficientsFromLibrary(const std::string &species);
 
    void run();
 
 
    void run(const std::vector<uint> &lUUIDs);
 
 
    EmpiricalModelCoefficients getEmpiricalModelCoefficients(uint UUID);
 
 
    FarquharModelCoefficients getFarquharModelCoefficients(uint UUID);
 
 
    C4ModelCoefficients getC4ModelCoefficients(uint UUID);
 
    void disableMessages();
 
    void enableMessages();
 
 
    void optionalOutputPrimitiveData(const char *label);
 
 
    void setCi(float Ci, const std::vector<uint> &UUIDs);
 
 
    void setCm(float Cm, const std::vector<uint> &UUIDs);
 
    void printDefaultValueReport() const;
 
 
    void printDefaultValueReport(const std::vector<uint> &UUIDs) const;
 
private:
    helios::Context *context;
 
    std::string model;
    EmpiricalModelCoefficients empiricalmodelcoeffs;
    FarquharModelCoefficients farquharmodelcoeffs;
    C4ModelCoefficients c4modelcoeffs;
 
    std::unordered_map<uint, EmpiricalModelCoefficients> empiricalmodel_coefficients;
    std::unordered_map<uint, FarquharModelCoefficients> farquharmodel_coefficients;
    std::unordered_map<uint, C4ModelCoefficients> c4model_coefficients;
 
    // Cache to avoid repeated Context lookups during run()
    mutable std::unordered_map<uint, EmpiricalModelCoefficients> material_coefficient_cache_empirical;
    mutable std::unordered_map<uint, FarquharModelCoefficients> material_coefficient_cache_farquhar;
    mutable std::unordered_map<uint, C4ModelCoefficients> material_coefficient_cache_c4;
 
    // Helper to retrieve coefficients with caching
    EmpiricalModelCoefficients getCoefficientsForPrimitive_Empirical(uint UUID) const;
    FarquharModelCoefficients getCoefficientsForPrimitive_Farquhar(uint UUID) const;
    C4ModelCoefficients getCoefficientsForPrimitive_C4(uint UUID) const;
 
    std::unordered_map<uint, float> previous_Ci;
 
    std::unordered_map<uint, float> manual_Ci;
 
    std::unordered_map<uint, float> manual_Cm;
 
    float evaluateEmpiricalModel(const EmpiricalModelCoefficients &params, float i_PAR, float TL, float CO2, float gM);
 
    float evaluateFarquharModel(const FarquharModelCoefficients &params, float i_PAR, float TL, float CO2, float gM, float &Ci, float &Gamma, int &limitation_state, int &TPU_flag, float &J_over_Jmax, helios::WarningAggregator &warnings);
 
    float evaluateC4Model(const C4ModelCoefficients &params, float i_PAR, float TL, float CO2, float gM, float &Ci, float &Cm, float &Vp, int &limitation_state, helios::WarningAggregator &warnings);
 
    float evaluateCi_Empirical(const EmpiricalModelCoefficients &params, float Ci, float CO2, float fL, float fT, float Rd, float gM) const;
 
    static float evaluateCi_Farquhar(float Ci, std::vector<float> &variables, const void *parameters);
 
    static float evaluateCi_C4(float Ci, std::vector<float> &variables, const void *parameters);
 
    static float evaluateCm_C4(float Cm, std::vector<float> &variables, const C4ModelCoefficients &params);
 
    static float respondToTemperature(const PhotosyntheticTemperatureResponseParameters *temperatureResponseParameters, float T);
 
    float i_PAR_default;
    float TL_default;
    float CO2_default;
    float gM_default;
    float gH_default;
 
    bool message_flag = true;
 
    std::vector<std::string> output_prim_data;
};
 
#endif