docs(foursail): add hotspot and polar BRDF tutorials

docs/make.jl

@@ -6,7 +6,8 @@ using Literate, UnitfulEquivalences, Distributions
 @info "Building docs with CanopyOptics source" path=pathof(CanopyOptics)
 
 function build()
-    tutorials = ["bilambertian.jl", "specular.jl", "foursail.jl", "dielectric.jl"]
+    tutorials = ["bilambertian.jl", "specular.jl", "foursail.jl",
+                 "foursail_hotspot.jl", "foursail_brdf.jl", "dielectric.jl"]
     tutorials_paths = [joinpath(@__DIR__, "src", "pages", "tutorials", tutorial) for tutorial in tutorials]
 
     for tutorial in tutorials_paths

docs/src/pages/tutorials/foursail_brdf.jl

@@ -0,0 +1,158 @@
+# # 4SAIL Polar BRDF and the Angular NDVI Effect
+
+# Where [the principal-plane tutorial](foursail_hotspot.md) sweeps a
+# single line through view-angle space, this one runs 4SAIL over the
+# full `(VZA, VAZ)` half-sphere in a single batched call and uses the
+# result to look at the *angular* dependence of NDVI.
+#
+# Two ideas are worth highlighting:
+#
+# * [`FourSAILGeometrySet`](@ref) lets you evaluate 4SAIL for an entire
+#   directional grid in one go — no Julia-level loop over view angles.
+# * NDVI is not a single number for a canopy: even a vegetation-only
+#   scene with a fixed leaf spectrum shows several percent of angular
+#   spread, with a peak on the hotspot side.
+using CanopyOptics
+using CairoMakie
+
+# ## Scene
+const sza_deg     = 30.0
+const lai         = 4.0
+const lambda_nm   = [685.0, 800.0]
+const soil_albedo = 0.0
+const hotspot_h   = 0.05
+
+# Dense (VZA, VAZ) grid. `FourSAILGeometrySet` consumes flattened
+# (VZA, VAZ) vectors; we reshape the results back to a 2D grid after.
+vzas_deg = collect(0.0:5.0:80.0)
+vazs_deg = collect(0.0:5.0:355.0)
+n_vza, n_vaz = length(vzas_deg), length(vazs_deg)
+
+vza_all = repeat(vzas_deg; inner = n_vaz)
+raa_all = repeat(vazs_deg; outer = n_vza)   # 4SAIL's raa = azimuth difference
+
+# ## Leaf optics from PROSPECT-PRO
+leaf = CanopyOptics.LeafProspectProProperties{Float64}(
+    N = 1.5, Ccab = 40.0, Ccar = 8.0,
+    Canth = 0.0, Cbrown = 0.0,
+    Cw = 0.012, Cm = 0.009,
+    Cprot = 0.0, Ccbc = 0.0,
+)
+
+opti = CanopyOptics.createLeafOpticalStruct(400.0:1.0:2500.0)
+T_grid, R_grid = CanopyOptics.prospect(leaf, opti)
+λ_grid   = [Float64(v.val) for v in opti.λ]
+grid_idx = [argmin(abs.(λ_grid .- λ)) for λ in lambda_nm]
+R_leaf   = R_grid[grid_idx]
+T_leaf   = T_grid[grid_idx]
+
+# ## One batched 4SAIL call for the entire half-sphere
+LD    = CanopyOptics.planophile_leaves2(Float64)
+geoms = FourSAILGeometrySet(LD;
+    sza_deg = sza_deg,
+    vza_deg = vza_all,
+    raa_deg = raa_all,
+    quadrature = CanopyOptics.CanopyQuadrature(n_leaf = 64, n_azimuth = 16),
+)
+
+result = foursail(R_leaf, T_leaf, soil_albedo, geoms, lai;
+                  hotspot = hotspot_h, mode = :full)
+
+# `result.rsot` is (n_wavelength × n_directions). Reshape each
+# wavelength row back to (n_vza × n_vaz):
+brf = Array{Float64}(undef, length(lambda_nm), n_vza, n_vaz)
+for iλ in eachindex(lambda_nm)
+    brf[iλ, :, :] .= permutedims(reshape(result.rsot[iλ, :], n_vaz, n_vza))
+end
+
+R_red = @view brf[1, :, :]
+R_nir = @view brf[2, :, :]
+NDVI  = (R_nir .- R_red) ./ (R_nir .+ R_red)
+
+println("685 nm BRF extrema: ", extrema(R_red))
+println("800 nm BRF extrema: ", extrema(R_nir))
+println("NDVI extrema:        ", extrema(NDVI))
+
+# Anchor the angular NDVI effect to nadir so we can see *how much* it
+# varies across the view hemisphere relative to the most common
+# space-based reference geometry.
+nadir_ndvi = NDVI[1, 1]
+ΔNDVI = NDVI .- nadir_ndvi
+println("NDVI at nadir: ", nadir_ndvi,
+        "   ΔNDVI extrema (relative to nadir): ", extrema(ΔNDVI))
+
+# ## Polar BRDF plots
+#
+# Closing the azimuth loop (`hcat(..., grid[:, 1])`) avoids the visible
+# seam at `vaz = 0/360°` in the polar surface plots. The red star marks
+# the hotspot direction at `(vza = sza, vaz = 0°)`, which in 4SAIL's
+# convention is the same azimuth as the sun.
+θ_rad = deg2rad.(vcat(vazs_deg, 360.0))
+close_az(grid) = hcat(grid, grid[:, 1])
+
+fig = Figure(size = (1500, 700))
+
+for (iλ, λ) in pairs(lambda_nm)
+    grid = brf[iλ, :, :]
+    ax = PolarAxis(fig[1, 2(iλ-1) + 1];
+        title    = "$(Int(round(λ))) nm BRF (rsot)",
+        theta_0  = π / 2,
+        direction = -1,
+        rticks   = 0:20:80,
+        rticklabelsize = 10,
+        thetaticks = (deg2rad.(0:45:315), string.(0:45:315) .* "°"))
+    hm = surface!(ax, θ_rad, vzas_deg, close_az(grid)';
+                  colormap = :viridis, shading = NoShading)
+    scatter!(ax, [deg2rad(0.0)], [sza_deg];
+             color = :red, markersize = 12, marker = :star5,
+             strokecolor = :white, strokewidth = 1.0)
+    Colorbar(fig[1, 2(iλ-1) + 2], hm; label = "BRF")
+end
+
+Δlim = maximum(abs, ΔNDVI)
+ax_ndvi = PolarAxis(fig[2, 1];
+    title = "NDVI",
+    theta_0 = π / 2, direction = -1,
+    rticks = 0:20:80, rticklabelsize = 10,
+    thetaticks = (deg2rad.(0:45:315), string.(0:45:315) .* "°"))
+hm_ndvi = surface!(ax_ndvi, θ_rad, vzas_deg, close_az(NDVI)';
+                   colormap = :viridis, shading = NoShading)
+Colorbar(fig[2, 2], hm_ndvi; label = "NDVI")
+
+ax_delta = PolarAxis(fig[2, 3];
+    title = "NDVI − NDVI(nadir)",
+    theta_0 = π / 2, direction = -1,
+    rticks = 0:20:80, rticklabelsize = 10,
+    thetaticks = (deg2rad.(0:45:315), string.(0:45:315) .* "°"))
+hm_delta = surface!(ax_delta, θ_rad, vzas_deg, close_az(ΔNDVI)';
+                    colormap = :balance, colorrange = (-Δlim, Δlim),
+                    shading = NoShading)
+Colorbar(fig[2, 4], hm_delta; label = "ΔNDVI")
+
+for ax in (ax_ndvi, ax_delta)
+    scatter!(ax, [deg2rad(0.0)], [sza_deg];
+             color = :red, markersize = 12, marker = :star5,
+             strokecolor = :white, strokewidth = 1.0)
+end
+
+Label(fig[0, 1:4],
+      "4SAIL polar BRDF and NDVI angular effect " *
+      "(SZA = $(Int(sza_deg))°, LAI = $(lai), planophile LAD, soil albedo = $(soil_albedo), h = $(hotspot_h))";
+      fontsize = 16)
+fig
+
+# ## Reading the figure
+#
+# * The 685 nm BRDF is dim and dominated by the hotspot peak: most of
+#   the red signal *is* the directly back-scattered branch, because
+#   chlorophyll absorbs almost everything that scatters into the diffuse
+#   field.
+# * The 800 nm BRDF is much brighter and noticeably flatter, with a
+#   broad bowl shape — multiple scattering smears the angular contrast
+#   out across the upper hemisphere.
+# * The two effects combine in the NDVI panel: NDVI dips on the hotspot
+#   side (where the red BRF is locally lifted) and rises on the
+#   forward-scatter side. Even with a fixed leaf spectrum and a
+#   horizontally homogeneous canopy this is a several-percent angular
+#   variation — large enough to matter for cross-instrument comparison
+#   if view geometries differ.

docs/src/pages/tutorials/foursail_hotspot.jl

@@ -0,0 +1,169 @@
+# # 4SAIL Hotspot in the Solar Principal Plane
+
+# This tutorial walks through a red/NIR principal-plane BRF diagnostic
+# with [`foursail`](@ref), using a realistic PROSPECT-PRO leaf and the
+# [`FourSAILGeometrySet`](@ref) batched entry point. It shows three
+# things in one figure:
+#
+# * the full BRF `rsot` with hotspot off (`h = 0`)
+# * the full BRF `rsot` with a Kuusk-style hotspot turned on
+# * the direct/direct branch `rsost` from the same hotspot solve
+#
+# The shape `rsot = rsost + rsodt` is what makes 4SAIL convenient as a
+# diagnostic tool — the embedded `rsost` is the path that carries the
+# hotspot, so you can read off the hotspot increment without rerunning
+# the model in a reduced mode.
+using CanopyOptics
+using CairoMakie
+
+# ## Geometry — signed VZA in the principal plane
+#
+# The "principal plane" is the vertical plane containing both the sun
+# and the observer. We sweep view zenith from −80° to +80° and split
+# the sign by relative azimuth:
+#
+# * `signed_vza < 0` → `raa = 0°`: backscatter side, same azimuth as
+#   the sun, the side that carries the hotspot.
+# * `signed_vza > 0` → `raa = 180°`: forward-scatter side.
+#
+# (4SAIL's `raa_deg` argument is the *azimuth difference* between sun
+# and view directions, so RAA = 0° puts the observer in the same
+# half-plane as the sun.)
+const sza_deg = 30.0
+const lai     = 4.0
+const hotspot_h = 0.05
+const lambda_nm = [685.0, 800.0]
+
+signed_vza = collect(-80.0:2.0:80.0)
+vza_abs    = abs.(signed_vza)
+raa_deg    = ifelse.(signed_vza .< 0, 0.0, 180.0)
+
+LD    = CanopyOptics.planophile_leaves2(Float64)
+geoms = FourSAILGeometrySet(LD;
+    sza_deg = sza_deg,
+    vza_deg = vza_abs,
+    raa_deg = raa_deg,
+    quadrature = CanopyOptics.CanopyQuadrature(n_leaf = 64, n_azimuth = 16),
+)
+
+# ## Leaf optics from PROSPECT-PRO
+#
+# Evaluate PROSPECT once on a 1 nm grid in the 400–2500 nm range, then
+# pull the two requested wavelengths out by nearest-neighbor index.
+leaf = CanopyOptics.LeafProspectProProperties{Float64}(
+    N = 1.5, Ccab = 40.0, Ccar = 8.0,
+    Canth = 0.0, Cbrown = 0.0,
+    Cw = 0.012, Cm = 0.009,
+    Cprot = 0.0, Ccbc = 0.0,
+)
+
+opti = CanopyOptics.createLeafOpticalStruct(400.0:1.0:2500.0)
+T_grid, R_grid = CanopyOptics.prospect(leaf, opti)
+λ_grid   = [Float64(v.val) for v in opti.λ]
+grid_idx = [argmin(abs.(λ_grid .- λ)) for λ in lambda_nm]
+R_leaf   = R_grid[grid_idx]
+T_leaf   = T_grid[grid_idx]
+
+println("PROSPECT leaf optics at the two probe wavelengths:")
+for (λ, ρ, τ) in zip(lambda_nm, R_leaf, T_leaf)
+    println("  λ = ", λ, " nm   ρ = ", round(ρ, digits = 4),
+            "   τ = ", round(τ, digits = 4))
+end
+
+# ## Two 4SAIL solves — hotspot off, hotspot on
+#
+# `foursail` returns a [`FourSAILResult`](@ref). For the batched
+# `FourSAILGeometrySet` path each field is a (n_wavelength × n_angle)
+# matrix.
+soil_albedo = 0.0
+
+sail_nohot = foursail(R_leaf, T_leaf, soil_albedo, geoms, lai;
+                      hotspot = 0.0, mode = :full)
+sail_hot   = foursail(R_leaf, T_leaf, soil_albedo, geoms, lai;
+                      hotspot = hotspot_h, mode = :full)
+
+size(sail_hot.rsot)
+
+# Quick numeric readout at three reference angles: the hotspot
+# (`signed_vza = -sza`), nadir, and the symmetric forward-scatter angle.
+hotspot_idx = findfirst(==(-sza_deg), signed_vza)
+nadir_idx   = findfirst(==(0.0), signed_vza)
+forward_idx = findfirst(==(sza_deg), signed_vza)
+
+for (iλ, λ) in pairs(lambda_nm)
+    println("λ = ", λ, " nm")
+    for (label, idx) in (("hotspot (-sza)", hotspot_idx),
+                         ("nadir",          nadir_idx),
+                         ("forward (+sza)", forward_idx))
+        idx === nothing && continue
+        println("  ", rpad(label, 16),
+                "  rsot(h=0) = ",      round(sail_nohot.rsot[iλ, idx], digits = 4),
+                "  rsot(h=$(hotspot_h)) = ", round(sail_hot.rsot[iλ, idx], digits = 4),
+                "  rsost = ",          round(sail_hot.rsost[iλ, idx], digits = 4))
+    end
+end
+
+# ## Plot — full BRF and the hotspot increment
+#
+# Top row: BRF curves themselves (with and without hotspot, plus the
+# direct/direct branch `rsost`). Bottom row: the hotspot increment
+# `rsot(h) − rsot(0)` is concentrated near `vza = −sza`, exactly where
+# the camera looks back along the sun direction.
+fig = Figure(size = (1100, 700))
+colors = (sail_nohot = :dodgerblue3,
+          sail_hot   = :seagreen4,
+          sail_ss    = :darkorange3,
+          increment  = :purple4)
+
+for (iλ, λ) in pairs(lambda_nm)
+    ax = Axis(fig[iλ, 1];
+        xlabel = iλ == length(lambda_nm) ?
+                 "signed VZA in solar principal plane (deg)" : "",
+        ylabel = "BRF",
+        title  = "$(Int(round(λ))) nm — principal-plane BRF")
+    lines!(ax, signed_vza, sail_nohot.rsot[iλ, :];
+           color = colors.sail_nohot, linewidth = 2.2, linestyle = :dash,
+           label = "rsot, h = 0")
+    lines!(ax, signed_vza, sail_hot.rsot[iλ, :];
+           color = colors.sail_hot,   linewidth = 2.5,
+           label = "rsot, h = $(hotspot_h)")
+    lines!(ax, signed_vza, sail_hot.rsost[iλ, :];
+           color = colors.sail_ss,    linewidth = 2.0, linestyle = :dot,
+           label = "rsost (direct/direct), h = $(hotspot_h)")
+    vlines!(ax, [-sza_deg]; color = :gray30, linestyle = :dashdot, linewidth = 1.2,
+            label = "hotspot direction")
+    vlines!(ax, [0.0];      color = :gray60, linestyle = :dot,    linewidth = 1.0)
+    axislegend(ax; position = :lt, framevisible = false)
+
+    axd = Axis(fig[iλ, 2];
+        xlabel = iλ == length(lambda_nm) ?
+                 "signed VZA in solar principal plane (deg)" : "",
+        ylabel = "Δ BRF",
+        title  = "$(Int(round(λ))) nm — hotspot increment")
+    lines!(axd, signed_vza, sail_hot.rsot[iλ, :] .- sail_nohot.rsot[iλ, :];
+           color = colors.increment, linewidth = 2.4,
+           label = "rsot(h) − rsot(0)")
+    hlines!(axd, [0.0]; color = :gray70, linewidth = 1.0)
+    vlines!(axd, [-sza_deg]; color = :gray30, linestyle = :dashdot, linewidth = 1.2)
+    axislegend(axd; position = :rt, framevisible = false)
+end
+
+Label(fig[0, 1:2],
+      "4SAIL principal-plane BRF and hotspot increment " *
+      "(SZA = $(Int(sza_deg))°, LAI = $(lai), planophile LAD, soil albedo = $(soil_albedo))";
+      fontsize = 16)
+fig
+
+# ## Reading the figure
+#
+# * Without the hotspot (`h = 0`) the BRF is smooth across the principal
+#   plane. The slight back–forward asymmetry is just the LAD-weighted
+#   single-scattering shape.
+# * Turning the hotspot on lifts the curve sharply near `vza = -sza`.
+#   The increment is essentially zero on the forward-scatter side, so
+#   the multiple-scattering remainder `rsodt` is hardly perturbed — the
+#   hotspot lives in the `rsost` branch.
+# * `rsost` itself sits well below the full BRF in the NIR because most
+#   of the 800 nm signal arrives through diffuse multiple scattering,
+#   not single scattering. In the red band leaves absorb most of what
+#   they scatter once, so `rsost` and the full `rsot` are much closer.