constant_gs_with_hydraulics.jl

using Photosynthesis
using Photosynthesis: step, Plant
using Plots
using Photosynthesis.PlantEcoPhys_jl: RHtoVPD

# Helper functions for environmental conditions
function f_tair(; t_mid=13, high=30, low=20)
    t = 0:23
    y = @. low + (high-low)*sin((t+12-t_mid)*pi/24)^2
    return y
end

function f_ppfd(; t_start=6, t_end=20, ppfd_max=1500)
    y = zeros(24)
    t = t_start:t_end
    y[t_start+1:t_end+1] = @. sin((t-t_start)*pi/(t_end-t_start))
    return y * ppfd_max
end

# Simulation parameters
n_days = 5
gs_constant = 0.15  # mol m-2 s-1

println("Running $(n_days)-day simulation with constant gs = $gs_constant mol m-2 s-1")
println("Includes photosynthesis + hydraulics model\n")

# Create plant with both photosynthesis and hydraulics
plant = Plant()

# Storage for results
n_hours = n_days * 24
time_hours = Float64[]
A_vec = Float64[]
E_vec = Float64[]
P_leaf_vec = Float64[]  # Xylem water potential at leaf
PLC_leaf_vec = Float64[]  # Percent loss of conductivity at leaf
soil_water_vec = Float64[]

# Run simulation
hour = 0
plant_alive = true

for day in 1:n_days
    # Daily environmental conditions
    Tair = f_tair(; high=35, low=20)
    PPFD = f_ppfd(; ppfd_max=1500)
    VPD = RHtoVPD.(50, Tair)

    for hour_of_day in 1:24
        external = ExternalPhotosynParams(
            Tair=Tair[hour_of_day],
            VPD=VPD[hour_of_day],
            PPFD=PPFD[hour_of_day]
        )

        # Run one hour with constant gs (includes photosynthesis + hydraulics)
        A, dead = step(plant, external, gs_constant)

        if dead
            println("Plant died at day $day, hour $hour_of_day")
            println("  Final soil water ratio: $(plant.hydraulic_state.soil_water_ratio)")
            println("  Final leaf PLC: $(plant.hydraulic_state.PLC[end])")
            global plant_alive = false
            break
        end

        # Store results
        push!(time_hours, hour)
        push!(A_vec, A)

        # Get transpiration from the last calculation
        E, _ = PhotosynEB(plant.leaf, external, gs_constant)
        push!(E_vec, E)

        # Hydraulic state variables
        push!(P_leaf_vec, plant.hydraulic_state.P[end])  # Pressure at top of plant
        push!(PLC_leaf_vec, plant.hydraulic_state.PLC[end])
        push!(soil_water_vec, plant.hydraulic_state.soil_water_ratio)

        global hour += 1
    end

    !plant_alive && break

    println("Day $day complete - Soil water: $(round(plant.hydraulic_state.soil_water_ratio, digits=3)), Leaf PLC: $(round(plant.hydraulic_state.PLC[end], digits=3))")
end

if plant_alive
    println("\nPlant survived all $n_days days!")
    println("Final soil water ratio: $(round(plant.hydraulic_state.soil_water_ratio, digits=3))")
    println("Final leaf PLC: $(round(plant.hydraulic_state.PLC[end], digits=3))")
end

# Calculate totals
total_A = sum(max.(0, A_vec))
total_E = sum(max.(0, E_vec))
println("\nTotal carbon gain: $(round(total_A, digits=1)) umol CO2 m-2")
println("Total water loss: $(round(total_E, digits=1)) mmol H2O m-2")

# Create plots
p1 = plot(time_hours, A_vec,
    label=false,
    xlabel="Time (hours)",
    ylabel="A (umol CO2 m-2 s-1)",
    title="Photosynthesis",
    linewidth=2
)

p2 = plot(time_hours, E_vec,
    label=false,
    xlabel="Time (hours)",
    ylabel="E (mmol H2O m-2 s-1)",
    title="Transpiration",
    linewidth=2,
    color=:blue
)

p3 = plot(time_hours, P_leaf_vec,
    label=false,
    xlabel="Time (hours)",
    ylabel="Xylem P (MPa)",
    title="Leaf Water Potential",
    linewidth=2,
    color=:red
)

p4 = plot(time_hours, PLC_leaf_vec,
    label=false,
    xlabel="Time (hours)",
    ylabel="PLC (fraction)",
    title="Leaf Cavitation (PLC)",
    linewidth=2,
    color=:orange
)
hline!(p4, [0.7], linestyle=:dash, color=:black, label="Death threshold")

p5 = plot(time_hours, soil_water_vec,
    label=false,
    xlabel="Time (hours)",
    ylabel="Soil water ratio",
    title="Soil Water Content",
    linewidth=2,
    color=:green
)

# Combine plots
plot(p1, p2, p3, p4, p5,
    layout=(3,2),
    size=(1200, 900),
    plot_title="Plant Simulation: Constant gs=$gs_constant (with Hydraulics)"
)

savefig("constant_gs_with_hydraulics.png")
println("\nPlot saved as 'constant_gs_with_hydraulics.png'")