selfTest.cpp

#include "PhotosynthesisModel.h"
 
#include <filesystem>
#include <fstream>
#include <sstream>
 
#define DOCTEST_CONFIG_IMPLEMENT
#include <doctest.h>
#include "doctest_utils.h"
 
using namespace std;
using namespace helios;
 
float err_tol = 1e-3;
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar Model Type Setting") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelType_Farquhar());
 
    FarquharModelCoefficients coeffs = photomodel.getFarquharModelCoefficients(UUID);
    DOCTEST_CHECK(coeffs.Vcmax == -1.0f); // Default uninitialized value
    DOCTEST_CHECK(coeffs.Jmax == -1.0f); // Default uninitialized value
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Light Response Curve - Farquhar Model") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    float Qin[9] = {0, 50, 100, 200, 400, 800, 1200, 1500, 2000};
    std::vector<float> AQ_expected{-2.37932, 8.29752, 12.5566, 16.2075, 16.7448, 16.7448, 16.7448, 16.7448, 16.7448};
 
    // Set Farquhar model coefficients
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(78.5f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.12f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
 
    for (int i = 0; i < 9; i++) {
        DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", Qin[i]));
        DOCTEST_CHECK_NOTHROW(photomodel.run());
 
        float A;
        DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
        DOCTEST_CHECK(A == doctest::Approx(AQ_expected[i]).epsilon(err_tol));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel CO2 Response Curve - Farquhar Model") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    float CO2[9] = {100, 200, 300, 400, 500, 600, 700, 800, 1000};
    std::vector<float> ACi_expected{1.72714, 7.32672, 12.4749, 17.199, 21.5281, 25.4923, 28.4271, 29.656, 31.3813};
 
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(78.5f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.12f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.0f)); // High light
 
    for (int i = 0; i < 9; i++) {
        DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", CO2[i]));
        DOCTEST_CHECK_NOTHROW(photomodel.run());
 
        float A;
        DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
        DOCTEST_CHECK(A == doctest::Approx(ACi_expected[i]).epsilon(err_tol));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Temperature Response Curve - Farquhar Model") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    float TL[7] = {270, 280, 290, 300, 310, 320, 330};
    std::vector<float> AT_expected{3.87863, 8.74928, 14.4075, 17.199, 16.0928, 11.6431, 4.00821};
 
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(78.5f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.12f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.0f)); // High light
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 400.0f)); // Normal CO2
 
    for (int i = 0; i < 7; i++) {
        DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", TL[i]));
        DOCTEST_CHECK_NOTHROW(photomodel.run());
 
        float A;
        DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
        DOCTEST_CHECK(A == doctest::Approx(AT_expected[i]).epsilon(err_tol));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar gm setter overloads") {
    FarquharModelCoefficients fq;
 
    // Default is +infinity (sentinel: Cc ≡ Ci, no mesophyll limitation).
    DOCTEST_CHECK(std::isinf(fq.getMesophyllConductance_gmTempResponse().value_at_25C));
 
    SUBCASE("Constant gm") {
        fq.setMesophyllConductance_gm(0.5f);
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(0.5f));
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().dHa == doctest::Approx(0.f));
    }
    SUBCASE("Arrhenius gm") {
        fq.setMesophyllConductance_gm(0.8f, 50.f);
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(0.8f));
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().dHa == doctest::Approx(50.f));
    }
    SUBCASE("Peaked Arrhenius gm") {
        fq.setMesophyllConductance_gm(0.4f, 45.f, 35.f);
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(0.4f));
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().Topt == doctest::Approx(273.15f + 35.f));
    }
    SUBCASE("Peaked Arrhenius gm with explicit dHd") {
        fq.setMesophyllConductance_gm(0.6f, 49.f, 40.f, 420.f);
        DOCTEST_CHECK(fq.getMesophyllConductance_gmTempResponse().dHd == doctest::Approx(420.f));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar default gm = infinity preserves legacy A") {
    // Reruns the three existing C3 reference sweeps (Light, CO2, Temperature) with
    // the default-infinite gm. The solver must short-circuit to the legacy Cc ≡ Ci
    // path and return values bit-for-bit identical (to within 1e-5) to the existing
    // hardcoded expected values. This guards the short-circuit in evaluateCi_Farquhar.
    const float strict_tol = 1e-5f;
 
    auto make_coeffs = []() {
        FarquharModelCoefficients fcoeffs;
        fcoeffs.setVcmax(78.5f, 65.33f);
        fcoeffs.setJmax(150.f, 43.54f);
        fcoeffs.setRd(2.12f, 46.39f);
        fcoeffs.setQuantumEfficiency_alpha(0.45f);
        return fcoeffs;
    };
 
    SUBCASE("Light response") {
        Context context_test;
        uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
        PhotosynthesisModel photomodel(&context_test);
        float Qin[9] = {0, 50, 100, 200, 400, 800, 1200, 1500, 2000};
        std::vector<float> AQ_expected{-2.37932f, 8.29752f, 12.5566f, 16.2075f, 16.7448f, 16.7448f, 16.7448f, 16.7448f, 16.7448f};
        FarquharModelCoefficients fcoeffs = make_coeffs();
        photomodel.setModelCoefficients(fcoeffs);
        for (int i = 0; i < 9; i++) {
            context_test.setPrimitiveData(UUID, "radiation_flux_PAR", Qin[i]);
            photomodel.run();
            float A;
            context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
            DOCTEST_CHECK(A == doctest::Approx(AQ_expected[i]).epsilon(strict_tol));
        }
    }
    SUBCASE("CO2 response") {
        Context context_test;
        uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
        PhotosynthesisModel photomodel(&context_test);
        float CO2[9] = {100, 200, 300, 400, 500, 600, 700, 800, 1000};
        std::vector<float> ACi_expected{1.72714f, 7.32672f, 12.4749f, 17.199f, 21.5281f, 25.4923f, 28.4271f, 29.656f, 31.3813f};
        FarquharModelCoefficients fcoeffs = make_coeffs();
        photomodel.setModelCoefficients(fcoeffs);
        context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.0f);
        for (int i = 0; i < 9; i++) {
            context_test.setPrimitiveData(UUID, "air_CO2", CO2[i]);
            photomodel.run();
            float A;
            context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
            DOCTEST_CHECK(A == doctest::Approx(ACi_expected[i]).epsilon(strict_tol));
        }
    }
    SUBCASE("Temperature response") {
        Context context_test;
        uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
        PhotosynthesisModel photomodel(&context_test);
        float TL[7] = {270, 280, 290, 300, 310, 320, 330};
        std::vector<float> AT_expected{3.87863f, 8.74928f, 14.4075f, 17.199f, 16.0928f, 11.6431f, 4.00821f};
        FarquharModelCoefficients fcoeffs = make_coeffs();
        photomodel.setModelCoefficients(fcoeffs);
        context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.0f);
        context_test.setPrimitiveData(UUID, "air_CO2", 400.0f);
        for (int i = 0; i < 7; i++) {
            context_test.setPrimitiveData(UUID, "temperature", TL[i]);
            photomodel.run();
            float A;
            context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
            DOCTEST_CHECK(A == doctest::Approx(AT_expected[i]).epsilon(strict_tol));
        }
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar very large gm reproduces legacy") {
    // Force the finite-gm quadratic path (gm = 1e9, below the 1e6 short-circuit threshold
    // is gm > 1e6; 1e9 keeps it finite with std::isfinite but above the threshold — so to
    // actually enter the quadratic branch we use gm slightly below threshold, 9e5). The
    // quadratic branch should agree with the legacy branch to within 1e-3.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photo_legacy(&context_test);
    PhotosynthesisModel photo_gm(&context_test);
 
    FarquharModelCoefficients legacy_coeffs;
    legacy_coeffs.setVcmax(80.f, 65.33f);
    legacy_coeffs.setJmax(150.f, 43.54f);
    legacy_coeffs.setRd(2.0f, 46.39f);
    legacy_coeffs.setQuantumEfficiency_alpha(0.45f);
 
    FarquharModelCoefficients gm_coeffs = legacy_coeffs;
    gm_coeffs.setMesophyllConductance_gm(9.0e5f); // just below 1e6 guard → enters quadratic branch
 
    photo_legacy.setModelCoefficients(legacy_coeffs);
    photo_gm.setModelCoefficients(gm_coeffs);
 
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f);
    float Ci_sweep[5] = {200.f, 400.f, 600.f, 800.f, 1000.f};
    for (float CO2: Ci_sweep) {
        context_test.setPrimitiveData(UUID, "air_CO2", CO2);
        photo_legacy.run();
        float A_legacy;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A_legacy);
        photo_gm.run();
        float A_gm;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A_gm);
        DOCTEST_CHECK(A_gm == doctest::Approx(A_legacy).epsilon(1e-3));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar gm finite regression") {
    // Regression anchors for the finite-gm quadratic branch with gm = 0.2 mol/m^2/s/bar.
    // Expected A values are produced by this implementation — not independent reference
    // data — and serve only to detect future unintended drift.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(80.f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
    fcoeffs.setMesophyllConductance_gm(0.2f);
 
    photomodel.setModelCoefficients(fcoeffs);
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f);
 
    float CO2_sweep[6] = {100.f, 200.f, 400.f, 600.f, 800.f, 1000.f};
    // Regression anchors: values produced by the current quadratic implementation.
    // A is strictly less than the legacy (infinite-gm) value at each point — mesophyll
    // diffusion limitation reduces net assimilation as expected.
    float A_expected[6] = {1.43302f, 5.90503f, 14.192f, 21.6148f, 27.8598f, 30.2768f};
    for (int i = 0; i < 6; ++i) {
        context_test.setPrimitiveData(UUID, "air_CO2", CO2_sweep[i]);
        photomodel.run();
        float A;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
        DOCTEST_CHECK(A == doctest::Approx(A_expected[i]).epsilon(1e-3));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Empirical Model Coefficients") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    EmpiricalModelCoefficients emp_coeffs;
    emp_coeffs.Asat = 20.0f;
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(emp_coeffs));
 
    EmpiricalModelCoefficients retrieved_coeffs = photomodel.getEmpiricalModelCoefficients(UUID);
    DOCTEST_CHECK(retrieved_coeffs.Asat == doctest::Approx(20.0f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Empirical Model Coefficients with UUIDs") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    EmpiricalModelCoefficients emp_coeffs;
    emp_coeffs.Asat = 20.0f;
    std::vector<uint> UUIDs = {UUID};
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(emp_coeffs, UUIDs));
 
    EmpiricalModelCoefficients retrieved_coeffs = photomodel.getEmpiricalModelCoefficients(UUID);
    DOCTEST_CHECK(retrieved_coeffs.Asat == doctest::Approx(20.0f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar Model Coefficients with UUIDs") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    FarquharModelCoefficients farq_coeffs;
    farq_coeffs.Vcmax = 90.0f;
    std::vector<uint> UUIDs = {UUID};
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(farq_coeffs, UUIDs));
 
    FarquharModelCoefficients retrieved_coeffs = photomodel.getFarquharModelCoefficients(UUID);
    DOCTEST_CHECK(retrieved_coeffs.Vcmax == doctest::Approx(90.0f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Invalid Input Handling") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Set up model with valid coefficients
    FarquharModelCoefficients fcoeffs;
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
 
    // Set invalid inputs
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", -50.0f)); // Should be clipped to 0
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", 150.0f)); // Should be replaced with 300K default
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", -10.0f)); // Should be clipped to 0
 
    // Disable messages to avoid warning output about invalid inputs
    photomodel.disableMessages();
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    photomodel.enableMessages();
 
    float A_invalid;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A_invalid));
 
    // Compute expected A using default values (PAR=0, TL=300K, CO2=390ppm)
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 0.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", 300.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 390.0f));
 
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    float A_expected;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A_expected));
 
    DOCTEST_CHECK(A_invalid == doctest::Approx(A_expected).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Optional Output Primitive Data") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("Ci"));
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("limitation_state"));
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("Gamma_CO2"));
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("electron_transport_ratio"));
 
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 400.0f));
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    float Ci, Gamma, J_over_Jmax;
    int limitation_state;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Ci", Ci));
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "limitation_state", limitation_state));
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Gamma_CO2", Gamma));
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "electron_transport_ratio", J_over_Jmax));
 
    DOCTEST_CHECK(Ci != 0.0f);
    DOCTEST_CHECK(Gamma != 0.0f);
    DOCTEST_CHECK(J_over_Jmax >= 0.0f);
    DOCTEST_CHECK(J_over_Jmax <= 1.0f);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Print Default Value Report") {
    Context context_test;
    PhotosynthesisModel photomodel(&context_test);
 
    // Disable messages to avoid verbose output
    photomodel.disableMessages();
 
    // Capture stdout from default value report
    capture_cout cout_buffer;
    DOCTEST_CHECK_NOTHROW(photomodel.printDefaultValueReport());
 
    // Verify the report was generated
    std::string report_output = cout_buffer.get_captured_output();
    DOCTEST_CHECK_MESSAGE(report_output.find("Photosynthesis Model Default Value Report") != std::string::npos, "Default value report should be generated");
 
    // Re-enable messages
    photomodel.enableMessages();
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Empirical Model Type Setting") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelType_Empirical());
 
    // After setting empirical model, should be able to get empirical coefficients
    EmpiricalModelCoefficients coeffs = photomodel.getEmpiricalModelCoefficients(UUID);
    DOCTEST_CHECK(coeffs.Asat > 0); // Should have positive default values
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Vector Coefficients with Size Mismatch") {
    Context context_test;
    uint UUID1 = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    uint UUID2 = context_test.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    std::vector<FarquharModelCoefficients> coeffs_vector(2);
    std::vector<uint> UUIDs = {UUID1}; // Only one UUID, but 2 coefficients
 
    // Capture the warning message from stderr
    capture_cerr cerr_buffer;
 
    // This should print a warning and return without setting coefficients
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(coeffs_vector, UUIDs));
 
    // Verify we captured the expected warning message
    std::string captured_warnings = cerr_buffer.get_captured_output();
    DOCTEST_CHECK_MESSAGE(captured_warnings.find("number of model coefficients (2) does not match number of UUIDs (1)") != std::string::npos, "Size mismatch should produce warning message");
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Vector Coefficients Matching Size") {
    Context context_test;
    uint UUID1 = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    uint UUID2 = context_test.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    std::vector<FarquharModelCoefficients> coeffs_vector(2);
    coeffs_vector[0].setVcmax(80.0f);
    coeffs_vector[1].setVcmax(90.0f);
    std::vector<uint> UUIDs = {UUID1, UUID2};
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(coeffs_vector, UUIDs));
 
    FarquharModelCoefficients retrieved1 = photomodel.getFarquharModelCoefficients(UUID1);
    FarquharModelCoefficients retrieved2 = photomodel.getFarquharModelCoefficients(UUID2);
    DOCTEST_CHECK(retrieved1.getVcmaxTempResponse().value_at_25C == doctest::Approx(80.0f).epsilon(err_tol));
    DOCTEST_CHECK(retrieved2.getVcmaxTempResponse().value_at_25C == doctest::Approx(90.0f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Library Species") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Disable messages to avoid verbose output about setting coefficients
    photomodel.disableMessages();
 
    // Test getting coefficients from library
    DOCTEST_CHECK_NOTHROW(photomodel.getFarquharCoefficientsFromLibrary("Almond"));
 
    // Test setting coefficients from library for all primitives
    DOCTEST_CHECK_NOTHROW(photomodel.setFarquharCoefficientsFromLibrary("Almond"));
 
    // Test setting coefficients from library for specific UUIDs
    std::vector<uint> UUIDs = {UUID};
    DOCTEST_CHECK_NOTHROW(photomodel.setFarquharCoefficientsFromLibrary("Apple", UUIDs));
 
    // Re-enable messages
    photomodel.enableMessages();
 
    // Verify coefficients were set correctly
    FarquharModelCoefficients almond_coeffs = photomodel.getFarquharModelCoefficients(UUID);
    DOCTEST_CHECK(almond_coeffs.getVcmaxTempResponse().value_at_25C == doctest::Approx(101.08f).epsilon(err_tol)); // Apple Vcmax
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Message Control") {
    Context context_test;
    PhotosynthesisModel photomodel(&context_test);
 
    // Test disabling and enabling messages
    DOCTEST_CHECK_NOTHROW(photomodel.disableMessages());
    DOCTEST_CHECK_NOTHROW(photomodel.enableMessages());
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Default Value Reports") {
    Context context_test;
    uint UUID1 = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    uint UUID2 = context_test.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Disable messages to avoid verbose output
    photomodel.disableMessages();
 
    // Capture stdout from default value report for all primitives
    capture_cout cout_buffer_all;
    DOCTEST_CHECK_NOTHROW(photomodel.printDefaultValueReport());
    std::string report_all = cout_buffer_all.get_captured_output();
    DOCTEST_CHECK_MESSAGE(report_all.find("Photosynthesis Model Default Value Report") != std::string::npos, "Default value report should be generated for all primitives");
 
    // Capture stdout from default value report for specific UUIDs
    capture_cout cout_buffer_subset;
    std::vector<uint> UUIDs = {UUID1, UUID2};
    DOCTEST_CHECK_NOTHROW(photomodel.printDefaultValueReport(UUIDs));
    std::string report_subset = cout_buffer_subset.get_captured_output();
    DOCTEST_CHECK_MESSAGE(report_subset.find("Photosynthesis Model Default Value Report") != std::string::npos, "Default value report should be generated for specific UUIDs");
 
    // Re-enable messages
    photomodel.enableMessages();
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Run with Specific UUIDs") {
    Context context_test;
    uint UUID1 = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    uint UUID2 = context_test.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Set up some basic data
    context_test.setPrimitiveData(UUID1, "radiation_flux_PAR", 500.0f);
    context_test.setPrimitiveData(UUID2, "radiation_flux_PAR", 600.0f);
 
    // Disable messages to avoid default value reports
    photomodel.disableMessages();
 
    // Test running model on specific UUIDs
    std::vector<uint> UUIDs = {UUID1};
    DOCTEST_CHECK_NOTHROW(photomodel.run(UUIDs));
 
    float A1, A2;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID1, "net_photosynthesis", A1));
 
    // UUID2 shouldn't have been processed by run(UUIDs), but run() without arguments processes all
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID2, "net_photosynthesis", A2));
 
    // Re-enable messages
    photomodel.enableMessages();
 
    DOCTEST_CHECK(A1 != 0.0f);
    DOCTEST_CHECK(A2 != 0.0f);
}
 
DOCTEST_TEST_CASE("PhotosyntheticTemperatureResponseParameters Constructors") {
    // Test default constructor
    PhotosyntheticTemperatureResponseParameters params_default;
    DOCTEST_CHECK(params_default.value_at_25C == doctest::Approx(100.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_default.dHa == doctest::Approx(60.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_default.dHd == doctest::Approx(600.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_default.Topt == doctest::Approx(10000.0f).epsilon(err_tol));
 
    // Test single parameter constructor with value only
    PhotosyntheticTemperatureResponseParameters params_single(75.0f);
    DOCTEST_CHECK(params_single.value_at_25C == doctest::Approx(75.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_single.dHa == doctest::Approx(0.0f).epsilon(err_tol));
 
    // Test constructor with value and negative dHa (should trigger else branch)
    PhotosyntheticTemperatureResponseParameters params_negative(80.0f, -10.0f);
    DOCTEST_CHECK(params_negative.value_at_25C == doctest::Approx(80.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_negative.dHa == doctest::Approx(-10.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_negative.dHd == doctest::Approx(600.0f).epsilon(err_tol)); // Should use default
 
    // Test constructor with value, dHa, and optimum temperature
    PhotosyntheticTemperatureResponseParameters params_three(90.0f, 50.0f, 35.0f);
    DOCTEST_CHECK(params_three.value_at_25C == doctest::Approx(90.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_three.dHa == doctest::Approx(50.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_three.dHd == doctest::Approx(500.0f).epsilon(err_tol)); // 10 * dHa
    DOCTEST_CHECK(params_three.Topt == doctest::Approx(273.15f + 35.0f).epsilon(err_tol));
 
    // Test constructor with value, negative dHa, and optimum temperature (triggers else branch)
    PhotosyntheticTemperatureResponseParameters params_three_neg(85.0f, -15.0f, 32.0f);
    DOCTEST_CHECK(params_three_neg.value_at_25C == doctest::Approx(85.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_three_neg.dHa == doctest::Approx(-15.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_three_neg.dHd == doctest::Approx(600.0f).epsilon(err_tol)); // Should use default 600
    DOCTEST_CHECK(params_three_neg.Topt == doctest::Approx(273.15f + 32.0f).epsilon(err_tol));
 
    // Test full constructor with all parameters
    PhotosyntheticTemperatureResponseParameters params_full(95.0f, 45.0f, 30.0f, 400.0f);
    DOCTEST_CHECK(params_full.value_at_25C == doctest::Approx(95.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_full.dHa == doctest::Approx(45.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_full.dHd == doctest::Approx(400.0f).epsilon(err_tol));
    DOCTEST_CHECK(params_full.Topt == doctest::Approx(273.15f + 30.0f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("FarquharModelCoefficients Temperature Response Methods") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    FarquharModelCoefficients coeffs;
 
    // Test all the temperature response setter methods and getters
 
    // Test Jmax with optimum temperature (3 parameter version)
    coeffs.setJmax(150.0f, 43.5f, 35.0f);
    PhotosyntheticTemperatureResponseParameters jmax_params = coeffs.getJmaxTempResponse();
    DOCTEST_CHECK(jmax_params.value_at_25C == doctest::Approx(150.0f).epsilon(err_tol));
    DOCTEST_CHECK(jmax_params.dHa == doctest::Approx(43.5f).epsilon(err_tol));
 
    // Test Jmax with optimum temperature and specified rate of decrease (4 parameter version)
    coeffs.setJmax(160.0f, 45.0f, 40.0f, 500.0f);
    jmax_params = coeffs.getJmaxTempResponse();
    DOCTEST_CHECK(jmax_params.dHd == doctest::Approx(500.0f).epsilon(err_tol));
 
    // Test Vcmax with optimum temperature (3 parameter version)
    coeffs.setVcmax(100.0f, 65.0f, 32.0f);
    PhotosyntheticTemperatureResponseParameters vcmax_params = coeffs.getVcmaxTempResponse();
    DOCTEST_CHECK(vcmax_params.value_at_25C == doctest::Approx(100.0f).epsilon(err_tol));
 
    // Test Vcmax with optimum temperature and specified rate of decrease (4 parameter version)
    coeffs.setVcmax(110.0f, 70.0f, 38.0f, 600.0f);
    vcmax_params = coeffs.getVcmaxTempResponse();
    DOCTEST_CHECK(vcmax_params.dHd == doctest::Approx(600.0f).epsilon(err_tol));
 
    // Test TPU methods
    coeffs.setTPU(8.0f);
    PhotosyntheticTemperatureResponseParameters tpu_params = coeffs.getTPUTempResponse();
    DOCTEST_CHECK(tpu_params.value_at_25C == doctest::Approx(8.0f).epsilon(err_tol));
    DOCTEST_CHECK(coeffs.TPU_flag == 1);
 
    coeffs.setTPU(9.0f, 25.0f);
    tpu_params = coeffs.getTPUTempResponse();
    DOCTEST_CHECK(tpu_params.dHa == doctest::Approx(25.0f).epsilon(err_tol));
 
    coeffs.setTPU(10.0f, 30.0f, 33.0f);
    tpu_params = coeffs.getTPUTempResponse();
    DOCTEST_CHECK(tpu_params.Topt == doctest::Approx(273.15f + 33.0f).epsilon(err_tol));
 
    coeffs.setTPU(11.0f, 35.0f, 36.0f, 450.0f);
    tpu_params = coeffs.getTPUTempResponse();
    DOCTEST_CHECK(tpu_params.dHd == doctest::Approx(450.0f).epsilon(err_tol));
 
    // Test Rd methods
    coeffs.setRd(2.5f, 46.0f, 34.0f);
    PhotosyntheticTemperatureResponseParameters rd_params = coeffs.getRdTempResponse();
    DOCTEST_CHECK(rd_params.value_at_25C == doctest::Approx(2.5f).epsilon(err_tol));
 
    coeffs.setRd(2.8f, 48.0f, 37.0f, 480.0f);
    rd_params = coeffs.getRdTempResponse();
    DOCTEST_CHECK(rd_params.dHd == doctest::Approx(480.0f).epsilon(err_tol));
 
    // Test quantum efficiency alpha methods
    coeffs.setQuantumEfficiency_alpha(0.4f, 20.0f);
    PhotosyntheticTemperatureResponseParameters alpha_params = coeffs.getQuantumEfficiencyTempResponse();
    DOCTEST_CHECK(alpha_params.dHa == doctest::Approx(20.0f).epsilon(err_tol));
 
    coeffs.setQuantumEfficiency_alpha(0.5f, 25.0f, 30.0f);
    alpha_params = coeffs.getQuantumEfficiencyTempResponse();
    DOCTEST_CHECK(alpha_params.Topt == doctest::Approx(273.15f + 30.0f).epsilon(err_tol));
 
    coeffs.setQuantumEfficiency_alpha(0.6f, 30.0f, 32.0f, 350.0f);
    alpha_params = coeffs.getQuantumEfficiencyTempResponse();
    DOCTEST_CHECK(alpha_params.dHd == doctest::Approx(350.0f).epsilon(err_tol));
 
    // Test light response curvature theta methods
    coeffs.setLightResponseCurvature_theta(0.1f, 15.0f);
    PhotosyntheticTemperatureResponseParameters theta_params = coeffs.getLightResponseCurvatureTempResponse();
    DOCTEST_CHECK(theta_params.dHa == doctest::Approx(15.0f).epsilon(err_tol));
 
    coeffs.setLightResponseCurvature_theta(0.2f, 20.0f, 35.0f);
    theta_params = coeffs.getLightResponseCurvatureTempResponse();
    DOCTEST_CHECK(theta_params.Topt == doctest::Approx(273.15f + 35.0f).epsilon(err_tol));
 
    coeffs.setLightResponseCurvature_theta(0.3f, 25.0f, 38.0f, 400.0f);
    theta_params = coeffs.getLightResponseCurvatureTempResponse();
    DOCTEST_CHECK(theta_params.dHd == doctest::Approx(400.0f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Edge Cases and Error Conditions") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Disable messages to avoid warning output
    photomodel.disableMessages();
 
    // Test with unknown optional output primitive data
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("unknown_primitive"));
 
    // Test running model with no coefficients set (should use defaults)
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    // Test getting coefficients for empirical model when none are set
    EmpiricalModelCoefficients emp_default = photomodel.getEmpiricalModelCoefficients(UUID);
    DOCTEST_CHECK(emp_default.Asat == doctest::Approx(18.18f).epsilon(err_tol)); // Default value
 
    // Test getting Farquhar coefficients when none are set
    FarquharModelCoefficients farq_default = photomodel.getFarquharModelCoefficients(UUID);
    DOCTEST_CHECK(farq_default.Vcmax == doctest::Approx(-1.0f).epsilon(err_tol)); // Uninitialized value
 
    // Test with extreme input values that trigger warnings - capture stderr
    // Re-enable messages to allow warnings to be produced
    photomodel.enableMessages();
 
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", -100.0f); // Negative PAR
    context_test.setPrimitiveData(UUID, "temperature", 150.0f); // Very low temperature
    context_test.setPrimitiveData(UUID, "air_CO2", -50.0f); // Negative CO2
    context_test.setPrimitiveData(UUID, "moisture_conductance", -0.1f); // Negative moisture conductance
    context_test.setPrimitiveData(UUID, "boundarylayer_conductance", -1.0f); // Negative boundary layer conductance
 
    // Capture warnings from extreme conditions (expected behavior)
    capture_cerr cerr_buffer;
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    // Verify we captured convergence warnings (expected for these extreme conditions)
    std::string captured_warnings = cerr_buffer.get_captured_output();
    DOCTEST_CHECK_MESSAGE(captured_warnings.find("Photosynthesis model failed to converge") != std::string::npos, "Extreme conditions should produce convergence warnings");
 
    // Verify the model still produces reasonable output despite bad inputs
    float A;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
    DOCTEST_CHECK(!std::isnan(A)); // Should not be NaN
}
 
DOCTEST_TEST_CASE("FarquharModelCoefficients TPU Testing") {
    Context context_test;
    PhotosynthesisModel photomodel(&context_test);
 
    // Test that TPU flag is properly set and cleared
    FarquharModelCoefficients coeffs;
    DOCTEST_CHECK(coeffs.TPU_flag == 0); // Should start as 0
 
    // Test all TPU setting methods set the flag to 1
    coeffs.setTPU(5.0f);
    DOCTEST_CHECK(coeffs.TPU_flag == 1);
 
    coeffs.setTPU(6.0f, 20.0f);
    DOCTEST_CHECK(coeffs.TPU_flag == 1);
 
    coeffs.setTPU(7.0f, 25.0f, 30.0f);
    DOCTEST_CHECK(coeffs.TPU_flag == 1);
 
    coeffs.setTPU(8.0f, 30.0f, 35.0f, 400.0f);
    DOCTEST_CHECK(coeffs.TPU_flag == 1);
}
 
DOCTEST_TEST_CASE("Farquhar TPU Singularity Fix") {
    // Regression test: The TPU limitation formula Wp = 3*TPU*Ci/(Ci - Gamma*) has a
    // singularity at Ci = Gamma*. When Ci is forced near Gamma*, the old formulation
    // produces non-physical assimilation values (huge or negative). The fix moves the
    // TPU comparison to the assimilation level where A_TPU = 3*TPU - Rd is a constant.
 
    Context context_test;
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    // Use Olive coefficients (species with TPU, known to trigger the bug)
    FarquharModelCoefficients coeffs = photomodel.getFarquharCoefficientsFromLibrary("Olive");
    DOCTEST_CHECK(coeffs.TPU_flag == 1);
 
    // Gamma* at 300K ≈ exp(19.02 - 37.83 / (0.0083144598 * 300)) ≈ 56.6 ppm
    // Test Ci values that bracket the singularity
    std::vector<float> Ci_values = {40.0f, 50.0f, 55.0f, 56.0f, 57.0f, 58.0f, 60.0f, 80.0f, 150.0f, 250.0f};
    std::vector<float> A_values;
 
    for (float Ci_val : Ci_values) {
        uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
 
        photomodel.setModelCoefficients(coeffs, std::vector<uint>{UUID});
        context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 1500.0f);  // Full sun
        context_test.setPrimitiveData(UUID, "temperature", 300.0f);          // ~27C
        context_test.setPrimitiveData(UUID, "air_CO2", 400.0f);
        context_test.setPrimitiveData(UUID, "moisture_conductance", 0.3f);
 
        // Force Ci to a specific value to directly probe the singularity region
        photomodel.setCi(Ci_val, std::vector<uint>{UUID});
 
        photomodel.run(std::vector<uint>{UUID});
 
        float A;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
        A_values.push_back(A);
    }
 
    // Gamma* at 300K ≈ 56.6 ppm. Below this, A must be negative (plant is respiring).
    // The buggy code produces A=25+ at Ci=40 due to the Wp singularity.
    for (size_t i = 0; i < A_values.size(); i++) {
        DOCTEST_CHECK_MESSAGE(std::isfinite(A_values[i]),
            "Non-finite A at Ci=" << Ci_values[i] << ": A=" << A_values[i]);
        if (Ci_values[i] < 55.0f) {
            DOCTEST_CHECK_MESSAGE(A_values[i] < 0.0f,
                "A should be negative below compensation point: Ci=" << Ci_values[i]
                << " A=" << A_values[i]);
        }
    }
 
    // A should be monotonically increasing with Ci
    for (size_t i = 1; i < A_values.size(); i++) {
        DOCTEST_CHECK_MESSAGE(A_values[i] >= A_values[i-1] - 0.5f,
            "Non-monotonic A: A[Ci=" << Ci_values[i] << "]=" << A_values[i]
            << " < A[Ci=" << Ci_values[i-1] << "]=" << A_values[i-1]);
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Complex Farquhar Testing") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Test complex Farquhar model with TPU limitation
    FarquharModelCoefficients coeffs;
    coeffs.setVcmax(100.0f, 65.0f);
    coeffs.setJmax(200.0f, 43.0f);
    coeffs.setRd(2.0f, 46.0f);
    coeffs.setQuantumEfficiency_alpha(0.4f);
    coeffs.setTPU(8.0f); // Enable TPU limitation
 
    photomodel.setModelCoefficients(coeffs);
 
    // Set conditions that might trigger different limitation states
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 1800.0f);
    context_test.setPrimitiveData(UUID, "temperature", 298.15f); // 25C
    context_test.setPrimitiveData(UUID, "air_CO2", 300.0f);
 
    // Enable optional outputs to trigger more code paths
    photomodel.optionalOutputPrimitiveData("Ci");
    photomodel.optionalOutputPrimitiveData("limitation_state");
    photomodel.optionalOutputPrimitiveData("Gamma_CO2");
 
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    // Verify outputs exist
    float A, Ci, Gamma;
    int limitation_state;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Ci", Ci));
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Gamma_CO2", Gamma));
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "limitation_state", limitation_state));
 
    // Verify reasonable values
    DOCTEST_CHECK(A != 0.0f);
    DOCTEST_CHECK(Ci > 0.0f);
    DOCTEST_CHECK(Gamma > 0.0f);
    DOCTEST_CHECK(limitation_state >= 0);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Empirical Model Edge Cases") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Test empirical model with extreme conditions
    EmpiricalModelCoefficients empirical_coeffs;
    empirical_coeffs.Asat = 25.0f;
    empirical_coeffs.theta = 70.0f;
 
    photomodel.setModelCoefficients(empirical_coeffs);
 
    // Test with very high light
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 3000.0f);
    context_test.setPrimitiveData(UUID, "temperature", 310.0f); // High temperature
    context_test.setPrimitiveData(UUID, "air_CO2", 800.0f); // High CO2
 
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    float A;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
    DOCTEST_CHECK(!std::isnan(A));
 
    // Test with zero light
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 0.0f);
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
    DOCTEST_CHECK(A < 0.0f); // Should be negative (respiration only)
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Temperature Response Edge Cases") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Test with temperature response that has dHa = 0 (constant response)
    FarquharModelCoefficients coeffs;
    coeffs.setVcmax(100.0f); // Single parameter - constant (dHa = 0)
    coeffs.setJmax(200.0f); // Single parameter - constant (dHa = 0)
    coeffs.setRd(2.0f); // Single parameter - constant (dHa = 0)
    coeffs.setQuantumEfficiency_alpha(0.4f); // Single parameter - constant (dHa = 0)
 
    photomodel.setModelCoefficients(coeffs);
 
    // Test at different temperatures - should give same result due to dHa = 0
    std::vector<float> temperatures = {280.0f, 300.0f, 320.0f};
    std::vector<float> results;
 
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 800.0f);
    context_test.setPrimitiveData(UUID, "air_CO2", 400.0f);
 
    for (float temp: temperatures) {
        context_test.setPrimitiveData(UUID, "temperature", temp);
        DOCTEST_CHECK_NOTHROW(photomodel.run());
 
        float A;
        DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "net_photosynthesis", A));
        results.push_back(A);
    }
 
    // Even with dHa = 0, there are other temperature effects (Kc, Ko, Gamma, etc.)
    // So just verify all computations completed without error and gave finite results
    for (float result: results) {
        DOCTEST_CHECK(!std::isnan(result));
        DOCTEST_CHECK(std::isfinite(result));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel - Material-Based Coefficients") {
    Context context;
    PhotosynthesisModel photomodel(&context);
 
    // Suppress messages
    photomodel.disableMessages();
 
    SUBCASE("Farquhar Model - Set and Retrieve Coefficients via Material") {
        // Create material
        context.addMaterial("test_leaf");
 
        // Set custom coefficients using setter methods
        FarquharModelCoefficients custom_coeffs;
        custom_coeffs.setVcmax(150.0f);
        custom_coeffs.setJmax(200.0f);
        custom_coeffs.setRd(2.0f);
        custom_coeffs.setQuantumEfficiency_alpha(0.4f);
        custom_coeffs.O = 210.0f;
 
        photomodel.setModelCoefficients("test_leaf", custom_coeffs);
 
        // Create primitives with this material
        uint p1 = context.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
        uint p2 = context.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
        context.assignMaterialToPrimitive(p1, "test_leaf");
        context.assignMaterialToPrimitive(p2, "test_leaf");
 
        // Verify material has the data
        DOCTEST_CHECK(context.doesMaterialDataExist("test_leaf", "photo_fq_Vcmax"));
        DOCTEST_CHECK(context.doesMaterialDataExist("test_leaf", "photo_fq_Jmax"));
        DOCTEST_CHECK(context.doesMaterialDataExist("test_leaf", "photo_fq_Rd"));
        DOCTEST_CHECK(context.doesMaterialDataExist("test_leaf", "photo_fq_alpha"));
 
        // Verify values
        float Vcmax, Jmax, Rd, alpha;
        context.getMaterialData("test_leaf", "photo_fq_Vcmax", Vcmax);
        context.getMaterialData("test_leaf", "photo_fq_Jmax", Jmax);
        context.getMaterialData("test_leaf", "photo_fq_Rd", Rd);
        context.getMaterialData("test_leaf", "photo_fq_alpha", alpha);
 
        DOCTEST_CHECK(Vcmax == doctest::Approx(150.0f));
        DOCTEST_CHECK(Jmax == doctest::Approx(200.0f));
        DOCTEST_CHECK(Rd == doctest::Approx(2.0f));
        DOCTEST_CHECK(alpha == doctest::Approx(0.4f));
    }
 
    SUBCASE("Multiple Primitives Share Material Coefficients") {
        // Create material
        context.addMaterial("shared_leaf");
 
        // Set coefficients using setter methods
        FarquharModelCoefficients coeffs;
        coeffs.setVcmax(120.0f);
        coeffs.setJmax(180.0f);
        coeffs.setRd(1.5f);
        coeffs.setQuantumEfficiency_alpha(0.35f);
        photomodel.setModelCoefficients("shared_leaf", coeffs);
 
        // Create many primitives with same material
        std::vector<uint> primitives;
        for (int i = 0; i < 50; i++) {
            uint p = context.addPatch(make_vec3(i, 0, 0), make_vec2(1, 1));
            primitives.push_back(p);
            context.assignMaterialToPrimitive(p, "shared_leaf");
        }
 
        // Verify all primitives have access to same coefficients via material
        uint shared_mat_id = context.getMaterialIDFromLabel("shared_leaf");
        for (uint p: primitives) {
            DOCTEST_CHECK(context.getPrimitiveMaterialID(p) == shared_mat_id);
        }
 
        // Material should only store coefficients once
        DOCTEST_CHECK(context.doesMaterialDataExist("shared_leaf", "photo_fq_Vcmax"));
        DOCTEST_CHECK(context.doesMaterialDataExist("shared_leaf", "photo_fq_Jmax"));
    }
 
    SUBCASE("Library Integration with Materials") {
        // Create materials for different species
        context.addMaterial("almond_leaf");
        context.addMaterial("grape_leaf");
 
        // Set coefficients from library
        photomodel.setFarquharCoefficientsFromLibrary("Almond", "almond_leaf");
        photomodel.setFarquharCoefficientsFromLibrary("Grape", "grape_leaf");
 
        // Verify both materials have coefficient data
        DOCTEST_CHECK(context.doesMaterialDataExist("almond_leaf", "photo_fq_Vcmax"));
        DOCTEST_CHECK(context.doesMaterialDataExist("grape_leaf", "photo_fq_Vcmax"));
 
        // Coefficients should be different for different species
        float almond_Vcmax, grape_Vcmax;
        context.getMaterialData("almond_leaf", "photo_fq_Vcmax", almond_Vcmax);
        context.getMaterialData("grape_leaf", "photo_fq_Vcmax", grape_Vcmax);
        DOCTEST_CHECK(almond_Vcmax != grape_Vcmax);
    }
 
    SUBCASE("Empirical Model - Material-Based API") {
        context.addMaterial("empirical_mat");
 
        EmpiricalModelCoefficients emp;
        emp.Asat = 25.0f;
        emp.theta = 70.0f;
        emp.Topt = 305.0f;
        emp.kC = 0.9f;
 
        DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients("empirical_mat", emp));
        DOCTEST_CHECK(context.doesMaterialDataExist("empirical_mat", "photo_emp_Asat"));
        DOCTEST_CHECK(context.doesMaterialDataExist("empirical_mat", "photo_emp_theta"));
        DOCTEST_CHECK(context.doesMaterialDataExist("empirical_mat", "photo_emp_Topt"));
 
        // Verify retrieved values
        float Asat_retrieved;
        context.getMaterialData("empirical_mat", "photo_emp_Asat", Asat_retrieved);
        DOCTEST_CHECK(Asat_retrieved == doctest::Approx(25.0f));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel - Manual Ci Setting") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    // Set up Farquhar model
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(78.5f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.12f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
 
    // Enable Ci output
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("Ci"));
 
    // Set environmental conditions
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 800.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", 298.15f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 400.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "moisture_conductance", 0.3f));
 
    // Set manual Ci
    float manual_Ci_value = 280.0f;
    std::vector<uint> UUIDs = {UUID};
    DOCTEST_CHECK_NOTHROW(photomodel.setCi(manual_Ci_value, UUIDs));
 
    // Run and verify Ci is manual value
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float Ci_retrieved;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Ci", Ci_retrieved));
    DOCTEST_CHECK(Ci_retrieved == doctest::Approx(manual_Ci_value).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel - Manual Ci Persistence") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(78.5f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.12f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("Ci"));
 
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 800.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", 298.15f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 400.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "moisture_conductance", 0.3f));
 
    // Set manual Ci
    float manual_Ci_value = 250.0f;
    std::vector<uint> UUIDs = {UUID};
    DOCTEST_CHECK_NOTHROW(photomodel.setCi(manual_Ci_value, UUIDs));
 
    // Run multiple times
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float Ci_run1;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Ci", Ci_run1));
 
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float Ci_run2;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Ci", Ci_run2));
 
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float Ci_run3;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID, "Ci", Ci_run3));
 
    // All runs should use same manual Ci
    DOCTEST_CHECK(Ci_run1 == doctest::Approx(manual_Ci_value).epsilon(err_tol));
    DOCTEST_CHECK(Ci_run2 == doctest::Approx(manual_Ci_value).epsilon(err_tol));
    DOCTEST_CHECK(Ci_run3 == doctest::Approx(manual_Ci_value).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel - Invalid Ci Error Handling") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    std::vector<uint> UUIDs = {UUID};
 
    // Test negative Ci
    capture_cerr cerr_capture;
    DOCTEST_CHECK_THROWS_AS(photomodel.setCi(-100.0f, UUIDs), std::runtime_error);
 
    // Test NaN
    DOCTEST_CHECK_THROWS_AS(photomodel.setCi(std::numeric_limits<float>::quiet_NaN(), UUIDs), std::runtime_error);
 
    // Test infinity
    DOCTEST_CHECK_THROWS_AS(photomodel.setCi(std::numeric_limits<float>::infinity(), UUIDs), std::runtime_error);
 
    // Test invalid UUID
    std::vector<uint> invalid_UUIDs = {99999};
    DOCTEST_CHECK_THROWS_AS(photomodel.setCi(300.0f, invalid_UUIDs), std::runtime_error);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel - Mixed Manual and Automatic Ci") {
    Context context_test;
    uint UUID1 = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    uint UUID2 = context_test.addPatch(make_vec3(2, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
 
    FarquharModelCoefficients fcoeffs;
    fcoeffs.setVcmax(78.5f, 65.33f);
    fcoeffs.setJmax(150.f, 43.54f);
    fcoeffs.setRd(2.12f, 46.39f);
    fcoeffs.setQuantumEfficiency_alpha(0.45f);
    DOCTEST_CHECK_NOTHROW(photomodel.setModelCoefficients(fcoeffs));
    DOCTEST_CHECK_NOTHROW(photomodel.optionalOutputPrimitiveData("Ci"));
 
    // Set same conditions for both primitives
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID1, "radiation_flux_PAR", 800.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID2, "radiation_flux_PAR", 800.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID1, "temperature", 298.15f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID2, "temperature", 298.15f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID1, "air_CO2", 400.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID2, "air_CO2", 400.0f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID1, "moisture_conductance", 0.3f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID2, "moisture_conductance", 0.3f));
 
    // Set manual Ci only for UUID1
    float manual_Ci_value = 250.0f;
    std::vector<uint> manual_UUIDs = {UUID1};
    DOCTEST_CHECK_NOTHROW(photomodel.setCi(manual_Ci_value, manual_UUIDs));
 
    // Run both
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    // UUID1 should have manual Ci
    float Ci1;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID1, "Ci", Ci1));
    DOCTEST_CHECK(Ci1 == doctest::Approx(manual_Ci_value).epsilon(err_tol));
 
    // UUID2 should have calculated Ci (different from manual)
    float Ci2;
    DOCTEST_CHECK_NOTHROW(context_test.getPrimitiveData(UUID2, "Ci", Ci2));
    DOCTEST_CHECK(Ci2 != doctest::Approx(manual_Ci_value).epsilon(err_tol));
    DOCTEST_CHECK(Ci2 > 0.0f);
    DOCTEST_CHECK(Ci2 < 400.0f);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel - Topt Parameter Validation") {
    Context context_test;
    PhotosynthesisModel photomodel(&context_test);
    FarquharModelCoefficients coeffs;
 
    SUBCASE("Negative Topt should throw error") {
        // Test negative Topt for Vcmax
        DOCTEST_CHECK_THROWS_WITH_AS(coeffs.setVcmax(100.0f, 60.0f, -10.0f),
                                     "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = -10.000000 C. Please check that temperature is provided in units of Celsius, not Kelvin.",
                                     std::runtime_error);
 
        // Test negative Topt for Jmax
        DOCTEST_CHECK_THROWS_WITH_AS(coeffs.setJmax(200.0f, 50.0f, -5.0f),
                                     "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = -5.000000 C. Please check that temperature is provided in units of Celsius, not Kelvin.",
                                     std::runtime_error);
 
        // Test negative Topt for Rd
        DOCTEST_CHECK_THROWS_WITH_AS(coeffs.setRd(2.0f, 40.0f, -15.0f),
                                     "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = -15.000000 C. Please check that temperature is provided in units of Celsius, not Kelvin.",
                                     std::runtime_error);
 
        // Test negative Topt for TPU
        DOCTEST_CHECK_THROWS_WITH_AS(coeffs.setTPU(5.0f, 50.0f, -20.0f),
                                     "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = -20.000000 C. Please check that temperature is provided in units of Celsius, not Kelvin.",
                                     std::runtime_error);
 
        // Test negative Topt for alpha
        DOCTEST_CHECK_THROWS_WITH_AS(coeffs.setQuantumEfficiency_alpha(0.5f, 30.0f, -25.0f),
                                     "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = -25.000000 C. Please check that temperature is provided in units of Celsius, not Kelvin.",
                                     std::runtime_error);
 
        // Test negative Topt for theta
        DOCTEST_CHECK_THROWS_WITH_AS(coeffs.setLightResponseCurvature_theta(0.7f, 35.0f, -30.0f),
                                     "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot be negative. Received Topt = -30.000000 C. Please check that temperature is provided in units of Celsius, not Kelvin.",
                                     std::runtime_error);
    }
 
    SUBCASE("Topt > 100C should throw error (likely Kelvin instead of Celsius)") {
        // Test Topt > 100°C for Vcmax (e.g., user provided 310 K instead of 36.85 C)
        DOCTEST_CHECK_THROWS_WITH_AS(
                coeffs.setVcmax(100.0f, 60.0f, 310.0f),
                "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = 310.000000 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).",
                std::runtime_error);
 
        // Test Topt > 100°C for Jmax
        DOCTEST_CHECK_THROWS_WITH_AS(
                coeffs.setJmax(200.0f, 50.0f, 305.0f),
                "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = 305.000000 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).",
                std::runtime_error);
 
        // Test Topt > 100°C for Rd
        DOCTEST_CHECK_THROWS_WITH_AS(
                coeffs.setRd(2.0f, 40.0f, 300.0f),
                "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = 300.000000 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).",
                std::runtime_error);
 
        // Test Topt > 100°C for TPU
        DOCTEST_CHECK_THROWS_WITH_AS(
                coeffs.setTPU(5.0f, 50.0f, 315.0f),
                "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = 315.000000 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).",
                std::runtime_error);
 
        // Test Topt > 100°C for alpha
        DOCTEST_CHECK_THROWS_WITH_AS(
                coeffs.setQuantumEfficiency_alpha(0.5f, 30.0f, 320.0f),
                "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = 320.000000 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).",
                std::runtime_error);
 
        // Test Topt > 100°C for theta
        DOCTEST_CHECK_THROWS_WITH_AS(
                coeffs.setLightResponseCurvature_theta(0.7f, 35.0f, 325.0f),
                "ERROR (PhotosyntheticTemperatureResponseParameters): Optimal temperature cannot exceed 100 C. Received Topt = 325.000000 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).",
                std::runtime_error);
    }
 
    SUBCASE("Valid Topt values should work correctly") {
        // Test valid Topt values (0-100°C range)
        DOCTEST_CHECK_NOTHROW(coeffs.setVcmax(100.0f, 60.0f, 35.0f));
        DOCTEST_CHECK_NOTHROW(coeffs.setJmax(200.0f, 50.0f, 40.0f));
        DOCTEST_CHECK_NOTHROW(coeffs.setRd(2.0f, 40.0f, 30.0f));
        DOCTEST_CHECK_NOTHROW(coeffs.setTPU(5.0f, 50.0f, 33.0f));
        DOCTEST_CHECK_NOTHROW(coeffs.setQuantumEfficiency_alpha(0.5f, 30.0f, 25.0f));
        DOCTEST_CHECK_NOTHROW(coeffs.setLightResponseCurvature_theta(0.7f, 35.0f, 28.0f));
 
        // Verify the Topt values were set correctly (converted to Kelvin internally)
        PhotosyntheticTemperatureResponseParameters vcmax_params = coeffs.getVcmaxTempResponse();
        DOCTEST_CHECK(vcmax_params.Topt == doctest::Approx(273.15f + 35.0f).epsilon(err_tol));
 
        PhotosyntheticTemperatureResponseParameters jmax_params = coeffs.getJmaxTempResponse();
        DOCTEST_CHECK(jmax_params.Topt == doctest::Approx(273.15f + 40.0f).epsilon(err_tol));
 
        PhotosyntheticTemperatureResponseParameters rd_params = coeffs.getRdTempResponse();
        DOCTEST_CHECK(rd_params.Topt == doctest::Approx(273.15f + 30.0f).epsilon(err_tol));
 
        PhotosyntheticTemperatureResponseParameters tpu_params = coeffs.getTPUTempResponse();
        DOCTEST_CHECK(tpu_params.Topt == doctest::Approx(273.15f + 33.0f).epsilon(err_tol));
 
        PhotosyntheticTemperatureResponseParameters alpha_params = coeffs.getQuantumEfficiencyTempResponse();
        DOCTEST_CHECK(alpha_params.Topt == doctest::Approx(273.15f + 25.0f).epsilon(err_tol));
 
        PhotosyntheticTemperatureResponseParameters theta_params = coeffs.getLightResponseCurvatureTempResponse();
        DOCTEST_CHECK(theta_params.Topt == doctest::Approx(273.15f + 28.0f).epsilon(err_tol));
    }
 
    SUBCASE("Edge cases: 0C and 100C should be valid") {
        // Test boundary values
        DOCTEST_CHECK_NOTHROW(coeffs.setVcmax(100.0f, 60.0f, 0.0f)); // 0°C should be valid
        DOCTEST_CHECK_NOTHROW(coeffs.setJmax(200.0f, 50.0f, 100.0f)); // 100°C should be valid
 
        PhotosyntheticTemperatureResponseParameters vcmax_params = coeffs.getVcmaxTempResponse();
        DOCTEST_CHECK(vcmax_params.Topt == doctest::Approx(273.15f).epsilon(err_tol));
 
        PhotosyntheticTemperatureResponseParameters jmax_params = coeffs.getJmaxTempResponse();
        DOCTEST_CHECK(jmax_params.Topt == doctest::Approx(373.15f).epsilon(err_tol));
    }
}
 
// -----------------------------------------------------------------------------
// C4 Model tests — von Caemmerer (2021) JXB 72:6003
// -----------------------------------------------------------------------------
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 Model Type Setting and Defaults") {
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    DOCTEST_CHECK_NOTHROW(photomodel.setModelType_C4());
 
    // Retrieve default C4 coefficients — paper Table 1 defaults for Setaria viridis
    C4ModelCoefficients c4 = photomodel.getC4ModelCoefficients(UUID);
 
    DOCTEST_CHECK(c4.Kc_25 == doctest::Approx(1210.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.Ko_25 == doctest::Approx(292000.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.Kp_25 == doctest::Approx(82.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.gamma_star_25 == doctest::Approx(3.81679e-4f).epsilon(err_tol));
    DOCTEST_CHECK(c4.Om_25 == doctest::Approx(210000.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.alpha_psII_fraction == doctest::Approx(0.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.x_etr_partition == doctest::Approx(0.4f).epsilon(err_tol));
    DOCTEST_CHECK(c4.Vpr == doctest::Approx(80.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.Rm_frac == doctest::Approx(0.5f).epsilon(err_tol));
    DOCTEST_CHECK(c4.fcyc == doctest::Approx(0.45f).epsilon(err_tol));
    DOCTEST_CHECK(c4.H_J == doctest::Approx(3.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.H_Jcyc == doctest::Approx(3.4f).epsilon(err_tol));
    DOCTEST_CHECK(c4.gbs == doctest::Approx(0.003f).epsilon(err_tol));
    DOCTEST_CHECK(c4.ao == doctest::Approx(0.047f).epsilon(err_tol));
 
    // Default temperature responses
    DOCTEST_CHECK(c4.getVpmaxTempResponse().value_at_25C == doctest::Approx(200.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.getVcmaxTempResponse().value_at_25C == doctest::Approx(40.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.getJmaxTempResponse().value_at_25C == doctest::Approx(247.69f).epsilon(err_tol));
    DOCTEST_CHECK(c4.getRdTempResponse().value_at_25C == doctest::Approx(1.f).epsilon(err_tol));
    DOCTEST_CHECK(c4.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(1.f).epsilon(err_tol));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 Setter Overloads - Temperature Response") {
    C4ModelCoefficients c4;
 
    SUBCASE("Constant Vpmax") {
        c4.setVpmax(150.f);
        DOCTEST_CHECK(c4.getVpmaxTempResponse().value_at_25C == doctest::Approx(150.f));
        DOCTEST_CHECK(c4.getVpmaxTempResponse().dHa == doctest::Approx(0.f));
    }
    SUBCASE("Arrhenius Vcmax") {
        c4.setVcmax(35.f, 70.f);
        DOCTEST_CHECK(c4.getVcmaxTempResponse().value_at_25C == doctest::Approx(35.f));
        DOCTEST_CHECK(c4.getVcmaxTempResponse().dHa == doctest::Approx(70.f));
    }
    SUBCASE("Peaked Arrhenius Jmax") {
        c4.setJmax(220.f, 45.f, 40.f);
        DOCTEST_CHECK(c4.getJmaxTempResponse().value_at_25C == doctest::Approx(220.f));
        DOCTEST_CHECK(c4.getJmaxTempResponse().Topt == doctest::Approx(273.15f + 40.f));
    }
    SUBCASE("Peaked Arrhenius Rd with explicit dHd") {
        c4.setRd(1.2f, 60.f, 42.f, 400.f);
        DOCTEST_CHECK(c4.getRdTempResponse().dHd == doctest::Approx(400.f));
    }
    SUBCASE("gm setters") {
        c4.setMesophyllConductance_gm(0.5f);
        DOCTEST_CHECK(c4.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(0.5f));
        c4.setMesophyllConductance_gm(0.8f, 50.f);
        DOCTEST_CHECK(c4.getMesophyllConductance_gmTempResponse().dHa == doctest::Approx(50.f));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - von Caemmerer 2021 spreadsheet verification at T=25C") {
    // Direct A vs Cm sweep at T=25 C, I=2000 umol/m²/s incident PAR. Reference values come from
    // the spreadsheet "C4__model_setaria__11-06-2021.xlsm" (sheet "A vs cm", column V "Amin").
    // Each test point supplies a Cm directly via setCm(), so the comparison tests the quadratic
    // math (Ac, Aj) without involving the Ci = Cm + A/gm coupling or the stomatal balance.
    //
    // This test explicitly pins the legacy vC2021 proton stoichiometry (fcyc=0.3, H_Jcyc=2)
    // because the reference spreadsheet was generated before the Woodford et al. (2025) update.
    // The new struct defaults (fcyc=0.45, H_Jcyc=3.4) are tested elsewhere. Keeping this test
    // ensures the generalized z = (H_J·(1−fcyc) + H_Jcyc·fcyc)/(h·(1−fcyc)) formula reduces
    // exactly to the vC2021 Eq. 31 form z = (3 − fcyc)/(h·(1−fcyc)) at H_J=3, H_Jcyc=2.
 
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    // Use paper defaults with three spreadsheet-compatibility adjustments:
    //  - Jmax(25 °C) set exactly to the Gaussian-evaluated value 247.6892 (peaked-Arrhenius fit matches only at 25 °C).
    //  - Vpr raised above Vpmax to disable the Eq. 19 cap Vp = min(Vp_MM, Vpr). The reference spreadsheet omits
    //    this cap, so our computed Ac at high Cm would diverge if we applied it. End-user defaults keep Vpr=80
    //    per paper Table 1 — the separate "Vpr cap activates" test verifies that behavior.
    //  - fcyc and H_Jcyc set to the vC2021 Eq. 31 values (not Woodford 2025 defaults) to match the 2021 spreadsheet.
    C4ModelCoefficients c4;
    c4.setJmax(247.6891682624015f);
    c4.Vpr = 1.0e4f;
    c4.fcyc = 0.3f;     // vC2021 Table 1 value; Woodford 2025 updates to 0.45.
    c4.H_Jcyc = 2.0f;   // vC2021 Eq. 31 implicit value (PGR5-only cyclic flow); Woodford 2025 updates to 3.4.
 
    photomodel.setModelCoefficients(c4);
 
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", 298.15f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 400.f));
 
    struct C4TestPoint {
        float Cm;
        float A_expected;
    };
    // Values from spreadsheet sheet "A vs cm", column F (Cm) and V (Amin) — paper Table 1 defaults
    const std::vector<C4TestPoint> points = {
            {0.1f, -0.567410f}, // boundary: Vp ≈ 0, A ≈ -Rd, exercises the discriminant clamp at the low edge
            {5.f,   8.695318f},  // Ac-limited (small Cm)
            {10.f,  16.159569f},
            {20.f,  25.628532f},
            {30.f,  30.140492f},
            {40.f,  32.364615f},
            {50.f,  33.604065f},
            {60.f,  34.375010f},
            {80.f,  34.835899f}, // Aj-limited from here
            {100.f, 34.862438f},
            {150.f, 34.928111f},
            {250.f, 35.056593f},
            {500.f, 35.361352f},
            {780.f, 35.675665f},
    };
 
    for (const auto &pt: points) {
        photomodel.setCm(pt.Cm, {UUID});
        DOCTEST_CHECK_NOTHROW(photomodel.run());
        float A;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
        // 0.5% tolerance — tight since the only approximations are float precision and the Jmax
        // peaked-Arrhenius-vs-Gaussian form (matched exactly at 25 °C).
        DOCTEST_CHECK_MESSAGE(A == doctest::Approx(pt.A_expected).epsilon(0.005),
                              "Cm=" << pt.Cm << " expected A=" << pt.A_expected << " got " << A);
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - Woodford 2025 default A vs Cm regression at T=25C") {
    // Companion to the vC2021 spreadsheet test. Same A vs Cm sweep at 25 °C with the same
    // Jmax/Vpr adjustments, but using the new struct defaults (fcyc=0.45, H_Jcyc=3.4, H_J=3)
    // from Woodford et al. (2025).
    //
    // NOTE: under these parameters (Vpmax=200, Vcmax=40, Jmax_25=247.69, Vpr=10⁴, fcyc=0.45)
    // Ac always limits at every tested Cm — Aj_plateau ≈ 52.3 stays well above
    // Ac_plateau ≈ 37.3, so this test exercises the Ac quadratic only. Cross-validated
    // exactly against the patched Woodford 2025 spreadsheet at 25 °C; for Aj-limited
    // validation at non-25 °C see "Patched Woodford spreadsheet anchors" below.
 
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    C4ModelCoefficients c4;
    c4.setJmax(247.6891682624015f);
    c4.Vpr = 1.0e4f;
    // fcyc, H_J, H_Jcyc use the Woodford 2025 struct defaults (0.45, 3, 3.4).
 
    photomodel.setModelCoefficients(c4);
 
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "temperature", 298.15f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 400.f));
 
    struct C4TestPoint {
        float Cm;
        float A_expected;
    };
    const std::vector<C4TestPoint> points = {
            // All points are Ac-limited under these parameters — Aj never wins at 25 °C.
            // Values computed from the patched Woodford 2025 spreadsheet (see
            // C4_spreadsheet_audit.md and tests/C4_test_anchors.csv at T_C=25.0).
            {0.1f, -0.567410f},
            {5.f,   8.695318f},
            {10.f,  16.159569f},
            {20.f,  25.628532f},
            {30.f,  30.140492f},
            {40.f,  32.364615f},
            {50.f,  33.604065f},
            {60.f,  34.375010f},
            {80.f,  35.2687f},
            {100.f, 35.7661f},
            {150.f, 36.3845f},
            {250.f, 36.8441f},
            {500.f, 37.1729f},
            {780.f, 37.2918f},
    };
 
    for (const auto &pt: points) {
        photomodel.setCm(pt.Cm, {UUID});
        DOCTEST_CHECK_NOTHROW(photomodel.run());
        float A;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
        DOCTEST_CHECK_MESSAGE(A == doctest::Approx(pt.A_expected).epsilon(0.005),
                              "Cm=" << pt.Cm << " expected A=" << pt.A_expected << " got " << A);
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - setCm round-trip: Ci = Cm + A/gm") {
    // When Cm is supplied directly, Ci is back-computed. Verify the relationship.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    photomodel.setModelType_C4();
    photomodel.optionalOutputPrimitiveData("Ci");
    photomodel.optionalOutputPrimitiveData("Cm");
 
    context_test.setPrimitiveData(UUID, "temperature", 298.15f);
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f);
    context_test.setPrimitiveData(UUID, "air_CO2", 400.f);
 
    C4ModelCoefficients defaults;
    const float gm_25 = defaults.getMesophyllConductance_gmTempResponse().value_at_25C;
 
    for (float Cm_in: {20.f, 60.f, 120.f}) {
        photomodel.setCm(Cm_in, {UUID});
        DOCTEST_CHECK_NOTHROW(photomodel.run());
        float A, Ci_out, Cm_out;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
        context_test.getPrimitiveData(UUID, "Ci", Ci_out);
        context_test.getPrimitiveData(UUID, "Cm", Cm_out);
        DOCTEST_CHECK(Cm_out == doctest::Approx(Cm_in).epsilon(1e-5));
        DOCTEST_CHECK(Ci_out == doctest::Approx(Cm_in + A / gm_25).epsilon(1e-3));
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - Temperature response wires through to A") {
    // Regression guard: catches the bug class "T-response getting silently disconnected from the C4 solver".
    // The setter tests only verify struct storage; the spreadsheet test runs only at 25 °C where the
    // Arrhenius factor degenerates to 1 regardless of dHa. This test runs at T = 308 K (35 °C),
    // where every default temperature-response actually contributes.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    photomodel.setModelType_C4();
 
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f);
    context_test.setPrimitiveData(UUID, "air_CO2", 400.f);
    photomodel.setCm(80.f, {UUID});
 
    // 25 °C baseline
    context_test.setPrimitiveData(UUID, "temperature", 298.15f);
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float A_25;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_25);
 
    // 35 °C — well below the Jmax peak (Topt = 43 °C) but with all rate constants substantially elevated
    context_test.setPrimitiveData(UUID, "temperature", 308.15f);
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float A_35;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_35);
 
    // Directional check: the bug we want to catch is the response being disconnected (A_25 == A_35).
    // C4 is near its temperature optimum at 35 °C, so A should rise meaningfully.
    DOCTEST_CHECK_MESSAGE(A_35 > A_25 * 1.1f,
                          "T-response not wired through: A(35°C)=" << A_35 << " not meaningfully > A(25°C)=" << A_25);
 
    // Tight regression anchors captured from the current implementation. Updating these is a deliberate
    // act and signals a behavioral change in either the temperature-response framework or the C4 solver.
    DOCTEST_CHECK(A_25 == doctest::Approx(34.012f).epsilon(0.005));
    DOCTEST_CHECK(A_35 == doctest::Approx(60.687f).epsilon(0.005));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - Vpr cap activates at high Cm") {
    // Per paper Eq. 19, Vp = min(Cm·Vpmax/(Cm+Kp), Vpr). At Cm above roughly 55 μbar (with Vpmax=200, Kp=82),
    // Vp_MM exceeds the default Vpr=80 and the cap must activate. Ensure the reported Vp never exceeds Vpr.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    photomodel.setModelType_C4();
    photomodel.optionalOutputPrimitiveData("Vp");
 
    context_test.setPrimitiveData(UUID, "temperature", 298.15f);
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f);
    context_test.setPrimitiveData(UUID, "air_CO2", 400.f);
    photomodel.setCi(600.f, {UUID});
 
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float Vp_out = 0.f;
    context_test.getPrimitiveData(UUID, "Vp", Vp_out);
    C4ModelCoefficients defaults;
    DOCTEST_CHECK(Vp_out <= defaults.Vpr + 1e-4f);
    DOCTEST_CHECK(Vp_out == doctest::Approx(defaults.Vpr).epsilon(1e-3)); // should hit the cap exactly at this Ci
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - Patched Woodford spreadsheet anchors") {
    // Validate A across the (T, Cm) grid against C4_test_anchors.csv, generated by
    // doc/generate_woodford_anchors.py from the patched Woodford 2025 spreadsheet
    // (γ* sign corrected, forward J(T,I) curvature chain installed). See
    // plugins/photosynthesis/doc/C4_spreadsheet_audit.md for the spreadsheet patches.
    //
    // To isolate validation of the Ac/Aj quadratics and the γ*(T) sign correction
    // from the separate Jmax-form difference (Helios uses peaked-Arrhenius;
    // spreadsheet uses Gaussian — they match at 25 °C and Topt=43 °C by construction
    // but diverge a few percent at intermediate temperatures), this test pins
    // Helios's Jmax to the Gaussian-derived value at each test temperature using
    // setJmax(value) which gives a constant (no-Arrhenius) T-response. Helios's
    // own Jmax temperature response is exercised separately by the "Temperature
    // response wires through to A" test above.
    //
    // The CSV transitions naturally between Ac- and Aj-limited regimes:
    //   - At T=25 °C, all anchors are Ac-limited (Ac plateau ~37, Aj plateau ~52).
    //   - At T=35 °C, Aj limits from Cm ~50 onward.
    //   - At T=40 °C, Aj limits from Cm ~40 onward.
    // Comparing Helios's net_photosynthesis (= min(Ac, Aj)) against Amin therefore
    // covers both quadratics across the temperature grid.
 
    auto gaussian_Jmax = [](float T_C) {
        const float Topt = 43.f;
        const float Omega = 26.f;
        const float Jmax_at_Topt = 400.f; // Helios C4 default
        const float dT = (T_C - Topt) / Omega;
        return Jmax_at_Topt * std::exp(-dT * dT);
    };
 
    const std::filesystem::path csv_path = helios::resolveFilePath("plugins/photosynthesis/tests/C4_test_anchors.csv");
    std::ifstream csv(csv_path);
    DOCTEST_REQUIRE_MESSAGE(csv.good(), "Failed to open " << csv_path.string());
 
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    DOCTEST_CHECK_NOTHROW(photomodel.setModelType_C4());
 
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f));
    DOCTEST_CHECK_NOTHROW(context_test.setPrimitiveData(UUID, "air_CO2", 400.f));
 
    std::string line;
    std::getline(csv, line); // header
    DOCTEST_REQUIRE(line.substr(0, 8) == "scenario");
 
    int n_points = 0;
    while (std::getline(csv, line)) {
        if (line.empty()) continue;
        std::stringstream ss(line);
        std::string scenario, T_str, Cm_str, Ac_str, Aj_str, Amin_str;
        std::getline(ss, scenario, ',');
        std::getline(ss, T_str, ',');
        std::getline(ss, Cm_str, ',');
        std::getline(ss, Ac_str, ',');
        std::getline(ss, Aj_str, ',');
        std::getline(ss, Amin_str, ',');
        const float T_C = std::stof(T_str);
        const float Cm = std::stof(Cm_str);
        const float Amin_expected = std::stof(Amin_str);
 
        C4ModelCoefficients c4;
        c4.Vpr = 1.0e4f;                       // disable Vpr cap (matches spreadsheet)
        c4.setJmax(gaussian_Jmax(T_C));        // dHa = 0 → Helios sees this exact Jmax at T_C
        photomodel.setModelCoefficients(c4);
 
        photomodel.setCm(Cm, {UUID});
        context_test.setPrimitiveData(UUID, "temperature", T_C + 273.15f);
        DOCTEST_CHECK_NOTHROW(photomodel.run());
        float A;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
        DOCTEST_CHECK_MESSAGE(A == doctest::Approx(Amin_expected).epsilon(0.005),
                              "T=" << T_C << "°C Cm=" << Cm
                              << " expected A=" << Amin_expected << " got A=" << A);
        ++n_points;
    }
    DOCTEST_CHECK_MESSAGE(n_points >= 60, "Expected ≥60 anchor rows in CSV; got " << n_points);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - Material-data caching round-trip") {
    // Verify that C4 coefficients assigned via material label are correctly serialized to
    // primitive-material data, retrieved on the next run(), and produce the same A as supplying
    // the same coefficients via the per-UUID setModelCoefficients overload.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    // Build a custom C4 coefficient set deliberately different from defaults
    C4ModelCoefficients custom;
    custom.setVcmax(60.f, 70.f);
    custom.setJmax(300.f);
    custom.Vpr = 1.0e4f; // disable cap so Vcmax differences register at Cm=50
    custom.alpha_psII_fraction = 0.05f;
 
    // Assign via material label
    const std::string mat = "test_c4_material";
    context_test.addMaterial(mat);
    context_test.assignMaterialToPrimitive(UUID, mat);
    photomodel.setModelCoefficients(mat, custom);
 
    context_test.setPrimitiveData(UUID, "temperature", 298.15f);
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 2000.f / 4.57f);
    context_test.setPrimitiveData(UUID, "air_CO2", 400.f);
    photomodel.setCm(50.f, {UUID});
    DOCTEST_CHECK_NOTHROW(photomodel.run());
    float A_material;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_material);
 
    // Re-create scenario applying the same coefficients via per-UUID setter on a fresh context
    Context ctx2;
    uint UUID2 = ctx2.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel pm2(&ctx2);
    pm2.disableMessages();
    pm2.setModelCoefficients(custom, {UUID2});
    ctx2.setPrimitiveData(UUID2, "temperature", 298.15f);
    ctx2.setPrimitiveData(UUID2, "radiation_flux_PAR", 2000.f / 4.57f);
    ctx2.setPrimitiveData(UUID2, "air_CO2", 400.f);
    pm2.setCm(50.f, {UUID2});
    DOCTEST_CHECK_NOTHROW(pm2.run());
    float A_uuid;
    ctx2.getPrimitiveData(UUID2, "net_photosynthesis", A_uuid);
 
    // Two paths must agree to high precision (only difference is the serialization round-trip)
    DOCTEST_CHECK(A_material == doctest::Approx(A_uuid).epsilon(1e-4));
 
    // Round-trip the coefficients struct itself and spot-check a few fields
    C4ModelCoefficients retrieved = photomodel.getC4ModelCoefficients(UUID);
    DOCTEST_CHECK(retrieved.getVcmaxTempResponse().value_at_25C == doctest::Approx(60.f));
    DOCTEST_CHECK(retrieved.getVcmaxTempResponse().dHa == doctest::Approx(70.f));
    DOCTEST_CHECK(retrieved.getJmaxTempResponse().value_at_25C == doctest::Approx(300.f));
    DOCTEST_CHECK(retrieved.Vpr == doctest::Approx(1.0e4f));
    DOCTEST_CHECK(retrieved.alpha_psII_fraction == doctest::Approx(0.05f));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 - Coefficient assignment per-UUID") {
    Context context_test;
    uint UUID1 = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    uint UUID2 = context_test.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    C4ModelCoefficients c4_hi;
    c4_hi.setVcmax(80.f);
    c4_hi.Vpr = 1.0e4f; // disable Vpr cap so Vcmax differences propagate to A in the Ac-limited regime
    photomodel.setModelCoefficients(c4_hi, {UUID1});
 
    C4ModelCoefficients c4_lo;
    c4_lo.setVcmax(20.f);
    c4_lo.Vpr = 1.0e4f;
    photomodel.setModelCoefficients(c4_lo, {UUID2});
 
    DOCTEST_CHECK(photomodel.getC4ModelCoefficients(UUID1).getVcmaxTempResponse().value_at_25C == doctest::Approx(80.f));
    DOCTEST_CHECK(photomodel.getC4ModelCoefficients(UUID2).getVcmaxTempResponse().value_at_25C == doctest::Approx(20.f));
 
    // Behavioral check: run the model on both UUIDs and confirm net_photosynthesis differs.
    // At Cm = 50 with Vcmax = 80 vs 20, the Ac-limited rate should clearly diverge.
    context_test.setPrimitiveData(UUID1, "temperature", 298.15f);
    context_test.setPrimitiveData(UUID2, "temperature", 298.15f);
    context_test.setPrimitiveData(UUID1, "radiation_flux_PAR", 2000.f / 4.57f);
    context_test.setPrimitiveData(UUID2, "radiation_flux_PAR", 2000.f / 4.57f);
    context_test.setPrimitiveData(UUID1, "air_CO2", 400.f);
    context_test.setPrimitiveData(UUID2, "air_CO2", 400.f);
    photomodel.setCm(50.f, {UUID1, UUID2});
    DOCTEST_CHECK_NOTHROW(photomodel.run());
 
    float A1, A2;
    context_test.getPrimitiveData(UUID1, "net_photosynthesis", A1);
    context_test.getPrimitiveData(UUID2, "net_photosynthesis", A2);
    DOCTEST_CHECK(A1 > A2 + 1.f); // higher Vcmax should yield meaningfully higher A
}
 
// -----------------------------------------------------------------------------
// C4 species parameter library tests
// -----------------------------------------------------------------------------
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 Library - Parameter switching") {
    // Verify each library entry populates the expected characteristic fields and that entries
    // are distinguishable from one another. Also exercise case-insensitive key matching.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    SUBCASE("SetariaViridis_vC2021") {
        DOCTEST_CHECK_NOTHROW(photomodel.setC4CoefficientsFromLibrary("SetariaViridis_vC2021"));
        C4ModelCoefficients c = photomodel.getC4ModelCoefficients(UUID);
        DOCTEST_CHECK(c.getVcmaxTempResponse().value_at_25C == doctest::Approx(40.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getVcmaxTempResponse().dHa == doctest::Approx(78.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getVpmaxTempResponse().value_at_25C == doctest::Approx(200.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getVpmaxTempResponse().dHa == doctest::Approx(50.1f).epsilon(err_tol));
        DOCTEST_CHECK(c.getRdTempResponse().value_at_25C == doctest::Approx(0.4f).epsilon(err_tol));
        DOCTEST_CHECK(c.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(1.0f).epsilon(err_tol));
        DOCTEST_CHECK(c.Kc_25 == doctest::Approx(1210.f).epsilon(err_tol));
        DOCTEST_CHECK(c.Kp_25 == doctest::Approx(82.f).epsilon(err_tol));
        DOCTEST_CHECK(c.dH_Kc == doctest::Approx(64.2f).epsilon(err_tol));
        DOCTEST_CHECK(c.fcyc == doctest::Approx(0.45f).epsilon(err_tol));
        DOCTEST_CHECK(c.H_Jcyc == doctest::Approx(3.4f).epsilon(err_tol));
    }
 
    SUBCASE("GenericC4_vC2000") {
        DOCTEST_CHECK_NOTHROW(photomodel.setC4CoefficientsFromLibrary("GenericC4_vC2000"));
        C4ModelCoefficients c = photomodel.getC4ModelCoefficients(UUID);
        DOCTEST_CHECK(c.getVcmaxTempResponse().value_at_25C == doctest::Approx(60.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getVpmaxTempResponse().value_at_25C == doctest::Approx(120.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getJmaxTempResponse().value_at_25C == doctest::Approx(400.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getRdTempResponse().value_at_25C == doctest::Approx(1.0f).epsilon(err_tol));
        // Q10 → Arrhenius: Ea ≈ 61.6 kJ/mol for Q10=2.3, 51.2 for Q10=2.0
        DOCTEST_CHECK(c.getVcmaxTempResponse().dHa == doctest::Approx(61.6f).epsilon(err_tol));
        DOCTEST_CHECK(c.getRdTempResponse().dHa == doctest::Approx(51.2f).epsilon(err_tol));
        DOCTEST_CHECK(c.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(1.0e4f).epsilon(err_tol));
        DOCTEST_CHECK(c.Kp_25 == doctest::Approx(80.f).epsilon(err_tol));
        DOCTEST_CHECK(c.fcyc == doctest::Approx(0.f).epsilon(err_tol));
    }
 
    SUBCASE("Maize_Massad2007") {
        DOCTEST_CHECK_NOTHROW(photomodel.setC4CoefficientsFromLibrary("Maize_Massad2007"));
        C4ModelCoefficients c = photomodel.getC4ModelCoefficients(UUID);
        DOCTEST_CHECK(c.getVcmaxTempResponse().value_at_25C == doctest::Approx(60.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getVcmaxTempResponse().dHa == doctest::Approx(67.3f).epsilon(err_tol));
        // Peaked-Arrhenius conversion of Massad's (Ea, Hd, dS) → Topt
        DOCTEST_CHECK(c.getVcmaxTempResponse().Topt == doctest::Approx(273.15f + 32.3f).epsilon(err_tol));
        DOCTEST_CHECK(c.getVpmaxTempResponse().Topt == doctest::Approx(273.15f + 43.1f).epsilon(err_tol));
        DOCTEST_CHECK(c.getJmaxTempResponse().Topt == doctest::Approx(273.15f + 31.5f).epsilon(err_tol));
        // Bernacchi 2001 C3-derived Kc/Ko — the critical internal-consistency constraint for Massad
        DOCTEST_CHECK(c.Kc_25 == doctest::Approx(650.f).epsilon(err_tol));
        DOCTEST_CHECK(c.Ko_25 == doctest::Approx(450000.f).epsilon(err_tol));
        DOCTEST_CHECK(c.dH_Kc == doctest::Approx(79.43f).epsilon(err_tol));
        DOCTEST_CHECK(c.Kp_25 == doctest::Approx(80.f).epsilon(err_tol));
        DOCTEST_CHECK(c.getRdTempResponse().value_at_25C == doctest::Approx(0.6f).epsilon(err_tol));
        DOCTEST_CHECK(c.getMesophyllConductance_gmTempResponse().value_at_25C == doctest::Approx(1.0e4f).epsilon(err_tol));
        DOCTEST_CHECK(c.fcyc == doctest::Approx(0.f).epsilon(err_tol));
    }
 
    SUBCASE("Case-insensitive key matching") {
        DOCTEST_CHECK_NOTHROW(photomodel.getC4CoefficientsFromLibrary("setariaviridis_vc2021"));
        DOCTEST_CHECK_NOTHROW(photomodel.getC4CoefficientsFromLibrary("MAIZE_MASSAD2007"));
        DOCTEST_CHECK_NOTHROW(photomodel.getC4CoefficientsFromLibrary("GenericC4"));
    }
 
    SUBCASE("Entries are distinguishable") {
        // Key library invariant: the three entries must differ in at least one characteristic
        // field each, otherwise picking between them has no effect on model behavior.
        C4ModelCoefficients setaria = photomodel.getC4CoefficientsFromLibrary("SetariaViridis_vC2021");
        C4ModelCoefficients generic = photomodel.getC4CoefficientsFromLibrary("GenericC4_vC2000");
        C4ModelCoefficients maize = photomodel.getC4CoefficientsFromLibrary("Maize_Massad2007");
 
        DOCTEST_CHECK(setaria.getVcmaxTempResponse().value_at_25C != generic.getVcmaxTempResponse().value_at_25C);
        DOCTEST_CHECK(setaria.Kc_25 != maize.Kc_25); // 1210 vs 650 — the critical Massad/Bernacchi distinction
        DOCTEST_CHECK(setaria.fcyc != generic.fcyc); // 0.45 vs 0
        // Simple Arrhenius (no peak) stores Topt as a large sentinel (≥10000 K by construction of
        // PhotosyntheticTemperatureResponseParameters); peaked Arrhenius stores an actual temperature in K.
        DOCTEST_CHECK(generic.getVpmaxTempResponse().Topt > 9000.f); // generic = simple Arrhenius, sentinel Topt
        DOCTEST_CHECK(maize.getVpmaxTempResponse().Topt < 400.f);    // maize = peaked Arrhenius, Topt ≈ 316 K
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 Library - Reasonable photosynthesis values") {
    // For each library entry, run the full model (with Ci→Cm coupling and stomatal balance) under
    // realistic midday conditions and confirm net photosynthesis falls within a reasonable C4
    // range, and that a low-light run produces smaller, still-positive A.
    const std::vector<std::string> species = {"SetariaViridis_vC2021", "GenericC4_vC2000", "Maize_Massad2007"};
 
    for (const auto &sp: species) {
        Context context_test;
        uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
        PhotosynthesisModel photomodel(&context_test);
        photomodel.disableMessages();
        photomodel.setC4CoefficientsFromLibrary(sp);
        photomodel.optionalOutputPrimitiveData("limitation_state");
 
        context_test.setPrimitiveData(UUID, "temperature", 298.15f);
        context_test.setPrimitiveData(UUID, "air_CO2", 400.f);
 
        // High-light midday
        context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 1800.f / 4.57f);
        DOCTEST_CHECK_NOTHROW(photomodel.run());
        float A_high = 0.f;
        int limitation_state = 0;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A_high);
        context_test.getPrimitiveData(UUID, "limitation_state", limitation_state);
        // Reasonable-range bounds are deliberately loose — the goal is to catch gross errors
        // (A near zero, or A in the hundreds from a unit mix-up), not to pin exact values.
        // Current implementation produces 34.6 (Setaria), 40.5 (Generic), 40.8 (Maize).
        bool A_high_in_range = (A_high >= 10.f) && (A_high <= 60.f);
        bool lim_valid = (limitation_state == 1) || (limitation_state == 2);
        DOCTEST_CHECK_MESSAGE(A_high_in_range, sp << ": A at high PAR out of reasonable C4 range (" << A_high << ")");
        DOCTEST_CHECK_MESSAGE(lim_valid, sp << ": unexpected limitation_state " << limitation_state);
 
        // Low-light: should be positive and lower than high-light
        context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 200.f / 4.57f);
        DOCTEST_CHECK_NOTHROW(photomodel.run());
        float A_low = 0.f;
        context_test.getPrimitiveData(UUID, "net_photosynthesis", A_low);
        bool A_low_valid = (A_low > 0.f) && (A_low < A_high);
        DOCTEST_CHECK_MESSAGE(A_low_valid, sp << ": low-light A=" << A_low << " (should be 0<A<" << A_high << ")");
    }
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 Library - Temperature discrimination") {
    // Setaria's Vcmax is simple Arrhenius (monotonic) while Maize's Vcmax is peaked at ~32 °C.
    // Confirm both give positive A at 25 °C and 35 °C, and catch the degenerate case where
    // switching species produces identical A (would indicate the library didn't actually switch).
    auto run_once = [](const std::string &sp, float T_K) {
        Context ctx;
        uint UUID = ctx.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
        PhotosynthesisModel pm(&ctx);
        pm.disableMessages();
        pm.setC4CoefficientsFromLibrary(sp);
        ctx.setPrimitiveData(UUID, "temperature", T_K);
        ctx.setPrimitiveData(UUID, "air_CO2", 400.f);
        ctx.setPrimitiveData(UUID, "radiation_flux_PAR", 1800.f / 4.57f);
        pm.run();
        float A = 0.f;
        ctx.getPrimitiveData(UUID, "net_photosynthesis", A);
        return A;
    };
 
    float setaria_25 = run_once("SetariaViridis_vC2021", 298.15f);
    float setaria_35 = run_once("SetariaViridis_vC2021", 308.15f);
    float maize_25 = run_once("Maize_Massad2007", 298.15f);
    float maize_35 = run_once("Maize_Massad2007", 308.15f);
 
    // All four scenarios should produce positive A
    DOCTEST_CHECK(setaria_25 > 0.f);
    DOCTEST_CHECK(setaria_35 > 0.f);
    DOCTEST_CHECK(maize_25 > 0.f);
    DOCTEST_CHECK(maize_35 > 0.f);
 
    // The two species must differ at both temperatures (otherwise the library didn't switch)
    DOCTEST_CHECK(std::abs(setaria_25 - maize_25) > 0.5f);
    DOCTEST_CHECK(std::abs(setaria_35 - maize_35) > 0.5f);
 
    // Setaria: monotonic Arrhenius Vcmax → A should rise meaningfully from 25→35 °C
    DOCTEST_CHECK(setaria_35 > setaria_25);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel C4 Library - Unknown species raises error") {
    Context context_test;
    context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    DOCTEST_CHECK_THROWS_AS(photomodel.getC4CoefficientsFromLibrary("NotARealSpecies"), std::runtime_error);
    DOCTEST_CHECK_THROWS_AS(photomodel.setC4CoefficientsFromLibrary(""), std::runtime_error);
    // Typo should fail fast rather than silently fall back
    DOCTEST_CHECK_THROWS_AS(photomodel.getC4CoefficientsFromLibrary("Maiz"), std::runtime_error);
}
 
 
// ---------------------------------------------------------------------------
// Regression tests for issues found in the deep-dive review.
// Each of these FAILS on the pre-fix code.
// ---------------------------------------------------------------------------
 
DOCTEST_TEST_CASE("PhotosynthesisModel Empirical temperature response is applied") {
    // REGRESSION: the documented f_T term (Tmin/Topt/q/Tref) was declared and serialized
    // but never used in the assimilation equation, so temperature entered only via Rd.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 200.f);
 
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    EmpiricalModelCoefficients emp;
    photomodel.setModelCoefficients(emp);
 
    // Evaluate at the optimum and well beyond it. f_T peaks at Topt, so A must decline.
    context_test.setPrimitiveData(UUID, "temperature", emp.Topt);
    photomodel.run();
    float A_at_opt;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_at_opt);
 
    context_test.setPrimitiveData(UUID, "temperature", emp.Topt + 12.f);
    photomodel.run();
    float A_hot;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_hot);
 
    // Assimilation at the optimum must exceed assimilation 12 K above it.
    DOCTEST_CHECK(A_at_opt > A_hot);
 
    // Changing Topt alone must change the answer at fixed leaf temperature.
    context_test.setPrimitiveData(UUID, "temperature", 310.f);
    photomodel.run();
    float A_default_topt;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_default_topt);
 
    EmpiricalModelCoefficients emp_shifted;
    emp_shifted.Topt = 280.f;
    photomodel.setModelCoefficients(emp_shifted);
    photomodel.run();
    float A_shifted_topt;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_shifted_topt);
 
    DOCTEST_CHECK(std::abs(A_default_topt - A_shifted_topt) > 1e-4f);
 
    // Below Tmin, assimilation is clamped off (gross assimilation zero, so A = -Rd < 0).
    EmpiricalModelCoefficients emp_cold;
    photomodel.setModelCoefficients(emp_cold);
    context_test.setPrimitiveData(UUID, "temperature", emp_cold.Tmin - 5.f);
    photomodel.run();
    float A_cold;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_cold);
    DOCTEST_CHECK(A_cold < 0.f);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Empirical invalid temperature coefficients throw below Tmin") {
    // REGRESSION: the Tref>Tmin / non-degenerate-denominator validation sat *after* the `TL <= Tmin`
    // early return, so an invalid coefficient set was only diagnosed on timesteps where the leaf
    // happened to be warmer than Tmin. On a cold timestep the model silently returned a well-formed
    // A = -Rd instead, hiding the misconfiguration for part of a diurnal/seasonal run. The coefficients
    // are a property of the set, not of the current temperature, so the check must fire either way.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
 
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    photomodel.setModelType_Empirical();
 
    // Tref below Tmin -- e.g. Tref and Tmin accidentally swapped, or one entered in Celsius.
    EmpiricalModelCoefficients bad_coeffs;
    bad_coeffs.Tmin = 290.f;
    bad_coeffs.Tref = 280.f;
    photomodel.setModelCoefficients(bad_coeffs);
 
    // Leaf at or below Tmin: the path that previously bypassed the check entirely.
    context_test.setPrimitiveData(UUID, "temperature", bad_coeffs.Tmin - 5.f);
    DOCTEST_CHECK_THROWS(photomodel.run());
 
    // Above Tmin the same invalid set must of course still throw.
    context_test.setPrimitiveData(UUID, "temperature", bad_coeffs.Tmin + 5.f);
    DOCTEST_CHECK_THROWS(photomodel.run());
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar material round-trip preserves TPU and theta") {
    // REGRESSION: setModelCoefficients(material,...) serialized TPU_flag but not TPU itself,
    // and never serialized theta. TPU silently reverted to the struct default (2.0) while
    // TPU_flag stayed 1, and the non-rectangular light response silently became rectangular.
    Context context_test;
    context_test.addMaterial("fq_roundtrip_mat");
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    context_test.assignMaterialToPrimitive(UUID, "fq_roundtrip_mat");
 
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
 
    FarquharModelCoefficients fq;
    fq.setVcmax(74.5f, 76.1f);
    fq.setJmax(180.2f, 23.0f);
    fq.setRd(1.3f, 46.39f);
    fq.setQuantumEfficiency_alpha(0.304f);
    fq.setLightResponseCurvature_theta(0.601f);
    fq.setTPU(7.7f, 24.0f);
 
    photomodel.setModelCoefficients("fq_roundtrip_mat", fq);
    FarquharModelCoefficients rt = photomodel.getFarquharModelCoefficients(UUID);
 
    DOCTEST_CHECK(rt.TPU_flag == 1);
    DOCTEST_CHECK(rt.getTPUTempResponse().value_at_25C == doctest::Approx(7.7f).epsilon(1e-4));
    DOCTEST_CHECK(rt.getTPUTempResponse().dHa == doctest::Approx(24.0f).epsilon(1e-4));
    DOCTEST_CHECK(rt.getLightResponseCurvatureTempResponse().value_at_25C == doctest::Approx(0.601f).epsilon(1e-4));
    DOCTEST_CHECK(rt.getQuantumEfficiencyTempResponse().value_at_25C == doctest::Approx(0.304f).epsilon(1e-4));
    DOCTEST_CHECK(rt.getVcmaxTempResponse().value_at_25C == doctest::Approx(74.5f).epsilon(1e-4));
 
    // A material with TPU never set must round-trip with TPU limitation still disabled.
    context_test.addMaterial("fq_notpu_mat");
    uint UUID2 = context_test.addPatch(make_vec3(1, 0, 0), make_vec2(1, 1));
    context_test.assignMaterialToPrimitive(UUID2, "fq_notpu_mat");
 
    FarquharModelCoefficients fq_notpu;
    fq_notpu.setVcmax(100.f, 65.33f);
    fq_notpu.setJmax(200.f, 46.36f);
    fq_notpu.setRd(1.5f, 46.39f);
    fq_notpu.setQuantumEfficiency_alpha(0.5f);
    photomodel.setModelCoefficients("fq_notpu_mat", fq_notpu);
 
    FarquharModelCoefficients rt2 = photomodel.getFarquharModelCoefficients(UUID2);
    DOCTEST_CHECK(rt2.TPU_flag == 0);
    DOCTEST_CHECK(rt2.getLightResponseCurvatureTempResponse().value_at_25C == doctest::Approx(0.f).epsilon(1e-6));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel peaked temperature response rejects dHd <= dHa") {
    // REGRESSION: respondToTemperature computed logf(dHd/dHa - 1), which is NaN when
    // dHd <= dHa. The NaN propagated into net_photosynthesis with no error and no warning
    // (the A == 0 convergence check is false for NaN).
    DOCTEST_CHECK_THROWS_AS(PhotosyntheticTemperatureResponseParameters(100.f, 65.f, 30.f, 50.f), std::runtime_error);
    DOCTEST_CHECK_THROWS_AS(PhotosyntheticTemperatureResponseParameters(100.f, 65.f, 30.f, 65.f), std::runtime_error);
 
    FarquharModelCoefficients fq;
    DOCTEST_CHECK_THROWS_AS(fq.setVcmax(100.f, 65.f, 30.f, 50.f), std::runtime_error);
    DOCTEST_CHECK_THROWS_AS(fq.setJmax(200.f, 46.f, 30.f, 40.f), std::runtime_error);
    DOCTEST_CHECK_THROWS_AS(fq.setRd(1.5f, 46.f, 30.f, 46.f), std::runtime_error);
    DOCTEST_CHECK_THROWS_AS(fq.setTPU(7.7f, 24.f, 30.f, 20.f), std::runtime_error);
 
    C4ModelCoefficients c4;
    DOCTEST_CHECK_THROWS_AS(c4.setVpmax(200.f, 50.f, 30.f, 40.f), std::runtime_error);
    DOCTEST_CHECK_THROWS_AS(c4.setJmax(247.f, 77.9f, 43.f, 70.f), std::runtime_error);
 
    // Valid peaked responses (dHd > dHa) must still be accepted and produce finite results.
    DOCTEST_CHECK_NOTHROW(fq.setVcmax(100.f, 65.f, 30.f, 200.f));
 
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 400.f);
    context_test.setPrimitiveData(UUID, "temperature", 300.f);
 
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    FarquharModelCoefficients good;
    good.setVcmax(100.f, 65.f, 30.f, 650.f);
    good.setJmax(200.f, 46.f);
    good.setRd(1.5f, 46.39f);
    good.setQuantumEfficiency_alpha(0.5f);
    photomodel.setModelCoefficients(good);
    photomodel.run();
    float A;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A);
    DOCTEST_CHECK(std::isfinite(A));
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar setters take precedence over legacy fields") {
    // REGRESSION: evaluateCi_Farquhar gated on the legacy public scalars (Vcmax/Jmax/Rd/alpha)
    // before consulting the temperature-response objects, so assigning the public field
    // silently discarded a prior setter call AND switched to a different T-response equation.
    Context context_test;
    uint UUID = context_test.addPatch(make_vec3(0, 0, 0), make_vec2(1, 1));
    context_test.setPrimitiveData(UUID, "radiation_flux_PAR", 200.f);
    context_test.setPrimitiveData(UUID, "temperature", 298.15f);
 
    PhotosynthesisModel photomodel(&context_test);
    photomodel.disableMessages();
    photomodel.setCi(250.f, {UUID});
 
    FarquharModelCoefficients baseline;
    baseline.setVcmax(100.f);
    baseline.setJmax(200.f);
    baseline.setRd(1.5f);
    baseline.setQuantumEfficiency_alpha(0.5f);
    photomodel.setModelCoefficients(baseline);
    photomodel.run();
    float A_setter;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_setter);
 
    // Calling the setter AFTER touching the legacy field must make the setter authoritative.
    FarquharModelCoefficients overridden;
    overridden.Vcmax = 50.f;
    overridden.setVcmax(100.f);
    overridden.setJmax(200.f);
    overridden.setRd(1.5f);
    overridden.setQuantumEfficiency_alpha(0.5f);
    photomodel.setModelCoefficients(overridden);
    photomodel.run();
    float A_overridden;
    context_test.getPrimitiveData(UUID, "net_photosynthesis", A_overridden);
 
    DOCTEST_CHECK(A_overridden == doctest::Approx(A_setter).epsilon(1e-4));
 
    // The setter must also clear the legacy sentinel so the two representations cannot diverge.
    DOCTEST_CHECK(overridden.Vcmax < 0.f);
}
 
DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar library values match documentation") {
    // REGRESSION: the documented species table diverged substantially from the shipped
    // library values (e.g. Almond doc Vcmax25=72.6 vs code 105.9).
    PhotosynthesisModel photomodel(nullptr);
    photomodel.disableMessages();
 
    struct LibExpect {
        const char *species;
        float Vcmax25;
        float Jmax25;
        float Rd25;
        float alpha;
    };
 
    // Values as published in plugins/photosynthesis/doc/Photosynthesis.dox
    const std::vector<LibExpect> expected{
            {"Almond", 72.6f, 144.2f, 0.2f, 0.094f},
            {"Walnut", 81.6f, 201.9f, 0.9f, 0.362f},
            {"Pistachio", 101.8f, 223.f, 1.5f, 0.216f},
            {"Apple", 101.08f, 167.03f, 3.00f, 0.432f},
            {"Cherry", 75.65f, 129.06f, 2.12f, 0.404f},
            {"Pear", 107.69f, 176.71f, 1.51f, 0.274f},
            {"Prune", 75.88f, 129.41f, 1.65f, 0.402f},
            {"Grape", 74.5f, 180.2f, 1.3f, 0.304f},
            {"Olive", 75.9f, 170.4f, 1.9f, 0.398f},
            {"Toyon", 52.8f, 142.4f, 0.8f, 0.29f},
            {"Elderberry", 37.7f, 149.7f, 1.3f, 0.202f},
            {"Maple", 96.4f, 168.f, 0.1f, 0.077f},
            {"Redbud", 68.5f, 132.4f, 0.8f, 0.41f},
    };
 
    for (const auto &e: expected) {
        FarquharModelCoefficients c = photomodel.getFarquharCoefficientsFromLibrary(e.species);
        DOCTEST_INFO("species = " << e.species);
        DOCTEST_CHECK(c.getVcmaxTempResponse().value_at_25C == doctest::Approx(e.Vcmax25).epsilon(1e-3));
        DOCTEST_CHECK(c.getJmaxTempResponse().value_at_25C == doctest::Approx(e.Jmax25).epsilon(1e-3));
        DOCTEST_CHECK(c.getRdTempResponse().value_at_25C == doctest::Approx(e.Rd25).epsilon(1e-3));
        DOCTEST_CHECK(c.getQuantumEfficiencyTempResponse().value_at_25C == doctest::Approx(e.alpha).epsilon(1e-3));
    }
 
    // Species documented with a TPU column must actually enable TPU limitation.
    for (const char *sp: {"Almond", "Walnut", "Pistachio"}) {
        FarquharModelCoefficients c = photomodel.getFarquharCoefficientsFromLibrary(sp);
        DOCTEST_INFO("species = " << sp);
        DOCTEST_CHECK(c.TPU_flag == 1);
    }
 
    // Species with no TPU column in the docs must leave TPU limitation off.
    for (const char *sp: {"Apple", "Cherry", "Pear", "Prune"}) {
        FarquharModelCoefficients c = photomodel.getFarquharCoefficientsFromLibrary(sp);
        DOCTEST_INFO("species = " << sp);
        DOCTEST_CHECK(c.TPU_flag == 0);
    }
 
    // The bare "Pistachio" key must resolve to the female cultivar, not the male one.
    FarquharModelCoefficients pist_default = photomodel.getFarquharCoefficientsFromLibrary("Pistachio");
    FarquharModelCoefficients pist_female = photomodel.getFarquharCoefficientsFromLibrary("PistachioFemale");
    FarquharModelCoefficients pist_male = photomodel.getFarquharCoefficientsFromLibrary("PistachioMale");
 
    DOCTEST_CHECK(pist_default.getVcmaxTempResponse().value_at_25C == doctest::Approx(pist_female.getVcmaxTempResponse().value_at_25C).epsilon(1e-4));
    DOCTEST_CHECK(pist_default.getJmaxTempResponse().value_at_25C == doctest::Approx(pist_female.getJmaxTempResponse().value_at_25C).epsilon(1e-4));
    DOCTEST_CHECK(pist_default.TPU_flag == pist_female.TPU_flag);
 
    // The male cultivar is a distinct, unpeaked fit with no TPU limitation.
    DOCTEST_CHECK(pist_male.getVcmaxTempResponse().value_at_25C == doctest::Approx(154.17f).epsilon(1e-3));
    DOCTEST_CHECK(pist_male.getJmaxTempResponse().value_at_25C == doctest::Approx(243.20f).epsilon(1e-3));
    DOCTEST_CHECK(pist_male.getRdTempResponse().value_at_25C == doctest::Approx(2.05f).epsilon(1e-3));
    DOCTEST_CHECK(pist_male.getQuantumEfficiencyTempResponse().value_at_25C == doctest::Approx(0.335f).epsilon(1e-3));
    DOCTEST_CHECK(pist_male.TPU_flag == 0);
    DOCTEST_CHECK(pist_male.getVcmaxTempResponse().value_at_25C != doctest::Approx(pist_female.getVcmaxTempResponse().value_at_25C).epsilon(1e-4));
}
 
int PhotosynthesisModel::selfTest(int argc, char **argv) {
    return helios::runDoctestWithValidation(argc, argv);
}