```
using VPL
using Distributions, Plots
# Data types
module TreeTypes
import VPL
# Meristem
struct Meristem <: VPL.Node end
# Bud
struct Bud <: VPL.Node end
# Node
struct Node <: VPL.Node end
# BudNode
struct BudNode <: VPL.Node end
# Internode (needs to be mutable to allow for changes over time)
Base.@kwdef mutable struct Internode <: VPL.Node
::Float64 = 0.10 # Internodes start at 10 cm
lengthend
# Leaf
Base.@kwdef struct Leaf <: VPL.Node
::Float64 = 0.20 # Leaves are 20 cm long
length::Float64 = 0.1 # Leaves are 10 cm wide
widthend
# Graph-level variables
Base.@kwdef struct treeparams
::Float64 = 0.1
growth::Float64 = 0.25
budbreak::Float64 = 140.0
phyllotaxis::Float64 = 30.0
leaf_angle::Float64 = 45.0
branch_angleend
end
import .TreeTypes
# Create geometry + color for the internodes
function VPL.feed!(turtle::Turtle, i::TreeTypes.Internode, vars)
# Rotate turtle around the head to implement elliptical phyllotaxis
rh!(turtle, vars.phyllotaxis)
HollowCylinder!(turtle, length = i.length, height = i.length/15, width = i.length/15,
= true, color = RGB(0.5,0.4,0.0))
move return nothing
end
# Create geometry + color for the leaves
function VPL.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.leaf_angle)
# Generate the leaf
Ellipse!(turtle, length = l.length, width = l.width, move = false,
= RGB(0.2,0.6,0.2))
color # Rotate turtle back to original direction
ra!(turtle, vars.leaf_angle)
return nothing
end
# Insertion angle for the bud nodes
function VPL.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.branch_angle)
end
# Rules
= Rule(TreeTypes.Meristem, rhs = mer -> TreeTypes.Node() +
meristem_rule Bud(), TreeTypes.Leaf()) +
(TreeTypes.Internode() + TreeTypes.Meristem())
TreeTypes.function prob_break(bud)
# We move to parent node in the branch where the bud was created
= parent(bud)
node # We count the number of internodes between node and the first Meristem
# moving down the graph
= hasDescendent(node, condition = n -> data(n) isa TreeTypes.Meristem)
check, steps = Int(ceil(steps/2)) # Because it will count both the nodes and the internodes
steps # Compute probability of bud break and determine whether it happens
if check
= min(1.0, steps*vars(bud).budbreak)
prob return rand() < prob
# If there is no meristem, an error happened since the model does not allow for this
else
error("No meristem found in branch")
end
end
= Rule(TreeTypes.Bud,
branch_rule = prob_break,
lhs = bud -> TreeTypes.BudNode() + TreeTypes.Internode() + TreeTypes.Meristem())
rhs
function create_tree(origin, growth, budbreak, orientation)
= T(origin) + RH(orientation) + TreeTypes.Internode() + TreeTypes.Meristem()
axiom = Graph(axiom = axiom, rules = (meristem_rule, branch_rule),
tree = TreeTypes.treeparams(growth = growth, budbreak = budbreak))
vars return tree
end
= Query(TreeTypes.Internode)
getInternode
function elongate!(tree, query)
for x in apply(tree, query)
= x.length*(1.0 + vars(tree).growth)
x.length end
end
function growth!(tree, query)
elongate!(tree, query)
rewrite!(tree)
end
function simulate(tree, query, nsteps)
= deepcopy(tree)
new_tree for i in 1:nsteps
growth!(new_tree, query)
end
return new_tree
end
= [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0]
origins = [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0]
orientations = rand(LogNormal(-2, 0.3), 10, 10)
growths = rand(Beta(2.0, 10), 10, 10)
budbreaks = vec(create_tree.(origins, growths, budbreaks, orientations));
forest = [simulate(tree, getInternode, 15) for tree in forest];
newforest render(newforest)
```

# Using ray tracer for ground cover

We use the previous simple model of the forest

To use the ray tracer to measure ground cover we just need to use `Black`

materials for the leaves, add a soil tile (also with a black material) and use a directional radiation source that is vertical. The ratio between absorbed power and emmitted power in the soil will give you the ground cover.

Since we are going to be using different types of ray tracers, we wil use the turtle message to use the right material in each case.

```
function choose_material(message)
if message == :original
nothing
elseif message == :ground_cover
Black()
else
error("The argument message must be :original or :ground_cover")
end
end
# Create geometry + color for the internodes
function VPL.feed!(turtle::Turtle, i::TreeTypes.Internode, vars)
# Rotate turtle around the head to implement elliptical phyllotaxis
rh!(turtle, vars.phyllotaxis)
HollowCylinder!(turtle, length = i.length, height = i.length/15, width = i.length/15,
= true, color = RGB(0.5,0.4,0.0),
move = choose_material(turtle.message))
material return nothing
end
# Create geometry + color for the leaves
function VPL.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.leaf_angle)
# Generate the leaf
Ellipse!(turtle, length = l.length, width = l.width, move = false,
= RGB(0.2,0.6,0.2), material = choose_material(turtle.message))
color # Rotate turtle back to original direction
ra!(turtle, vars.leaf_angle)
return nothing
end
```

We create a scene with the right message and add a soil tile:

```
= Scene(newforest, message = :ground_cover)
scene = Rectangle(length = 21.0, width = 21.0)
soil rotatey!(soil, pi/2)
translate!(soil, Vec(0.0, 10.5, 0.0))
VPL.= Black()
soil_mat add!(scene, mesh = soil, material = soil_mat) VPL.
```

We now add a single vertical source on top of the scene

`= DirectionalSource(scene, θ = 0.0, Φ = 0.0, radiosity = 1.0, nrays = 1_000_000) source `

And we run the ray tracer without grid cloner (no need for vertical source and black materials, just ignore the warning!)

```
= RTSettings(nx = 0, ny = 0, parallel = true)
settings = RayTracer(scene, source, settings = settings, acceleration = BVH,
rtobj = SAH{3}(5, 10));
rule trace!(rtobj)
```

And the ground cover is just the total power absorbed by the soil tile scaled by the emmitted power:

`1 - power(soil_mat)[1]/1_000_000/source.power[1]`