Bi-Lambertian Canopy Scattering

CanopyOptics represents the diffuse part of a leaf as a bi-Lambertian scatterer with leaf reflectance R and transmittance T.

using CanopyOptics
using CairoMakie
using Distributions
using Base64

Use a positive-μ quadrature grid and a standard planophile leaf-angle distribution.

μ, w = CanopyOptics.gauleg(12, 0.0, 1.0)
LD = CanopyOptics.planophile_leaves2()

LeafDistribution{Float64}(Beta{Float64}(α=1.0, β=2.768), 0.6366197723675814)

Build the leaf scattering model. The single-scattering albedo is ϖ = R + T; CanopyOptics divides that factor out of the returned Z matrices so downstream layer solvers can apply ϖ explicitly.

leaf = CanopyOptics.BiLambertianCanopyScattering(R = 0.4, T = 0.2)
quadrature = CanopyOptics.CanopyQuadrature(n_leaf = 64)

CanopyQuadrature(64, 48)

Closed-form Fourier stack

Passing a moment range to compute_Z_matrices returns all cosine Fourier moments in one call. The array layout is Z[i_out, j_in, m+1] for the common 0:m_max range.

m_max = 8
Z⁺⁺, Z⁻⁺ = CanopyOptics.compute_Z_matrices(leaf, μ, LD, 0:m_max; quadrature)
size(Z⁺⁺)

(12, 12, 9)

Z⁺⁺ is same-sign transmission and Z⁻⁺ is sign-change reflection. For conservative leaves, the m = 0 quadrature-weighted column integral of Z⁺⁺ + Z⁻⁺ is approximately 2 in this normalization.

column_flux = vec(sum(w .* (Z⁺⁺[:, :, 1] .+ Z⁻⁺[:, :, 1]), dims = 1))
extrema(column_flux)

(1.9999276998633733, 2.0001333845187514)

Single-moment calls slice the same analytic closure.

Z⁺⁺₂, Z⁻⁺₂ = CanopyOptics.compute_Z_matrices(leaf, μ, LD, 2; quadrature)
maximum(abs.(Z⁺⁺₂ .- Z⁺⁺[:, :, 3]))

0.0

A single m = 0 call slices the same analytic stack.

Z⁺⁺₀, Z⁻⁺₀ = CanopyOptics.compute_Z_matrices(leaf, μ, LD, 0; quadrature)
maximum(abs.(Z⁻⁺₀ .- Z⁻⁺[:, :, 1]))

0.0

Changing the leaf-angle distribution

Leaf-angle distributions can be supplied directly as distributions on [0, 1], with the standard 2 / π scaling from inclination angle to the normalized interval.

erectophile = CanopyOptics.LeafDistribution(Beta(5, 2), 2 / π)
Z_erect⁺⁺, Z_erect⁻⁺ =
    CanopyOptics.compute_Z_matrices(leaf, μ, erectophile, 0:m_max; quadrature)

maximum(abs.(Z_erect⁻⁺[:, :, 1] .- Z⁻⁺[:, :, 1]))

1.3205921849330804

Lightweight animation

CairoMakie renders without a Python backend. This animation sweeps a family of beta leaf-angle distributions and updates the m = 0 reflection matrix.

function reflection_matrix_for_beta(a, b)
    LD = CanopyOptics.LeafDistribution(Beta(a, b), 2 / π)
    _, Z⁻⁺ = CanopyOptics.compute_Z_matrices(leaf, μ, LD, 0:0; quadrature)
    return Z⁻⁺[:, :, 1]
end

αβ = [(0.8, 5.0), (1.5, 4.0), (2.5, 2.5), (4.0, 1.5), (5.0, 0.8)]
Zobs = Observable(reflection_matrix_for_beta(αβ[1]...))

fig = Figure(size = (560, 420))
ax = Axis(fig[1, 1]; xlabel = "outgoing μ", ylabel = "incoming μ",
          title = "Z⁻⁺, m = 0")
hm = heatmap!(ax, μ, μ, Zobs; colormap = :viridis)
Colorbar(fig[1, 2], hm; label = "Z⁻⁺")

path = record(fig, "anim_bilambertian.gif", eachindex(αβ); framerate = 3) do i
    Zobs[] = reflection_matrix_for_beta(αβ[i]...)
    ax.title = "Z⁻⁺, m = 0, Beta$(αβ[i])"
end
HTML("<img src=\"data:image/gif;base64,$(base64encode(read(path)))\" />")

This page was generated using Literate.jl.