This document shows, with all the gory details, how the calculations are made for finding the trees that intersect a given GLI light ray. Recall that light rays are defined by an origin point (x, y, z), an azimuth angle, and an altitude angle.
The code begins by calculating a maximum search radius for a given light ray. It takes the tallest possible tree (the largest value in the allometry file for maximum tree height) and then determines the distance at which the tree top just grazes the azimuth angle. (See diagrams below for how that relationship would be worked out.)
The code searches within this maximum radius for trees taller than the z coordinate of the light ray origin. All trees found in this way will be assessed as possible interceptors. (It would be more efficient, of course, to narrow the search based on azimuth direction. However, it's much easier to rip off SORTIE's circle-searching code.)
The code assesses azimuth intersection and altitude intersection separately.
When determining the azimuth, the code takes a two-dimensional projection of the tree and its neighbor on the forest floor.
In the figure above, the target tree is at (x1, y1) and the neighbor tree is at (x2, y2). The neighbor's canopy radius is r. θ is the azimuth angle of the light ray. To determine whether or not the circle of the neighbor’s canopy intersects the azimuth line, the code solves the circle-and-line system of two equations. The equation for a line is
To simplify the problem, the code puts the two trees in their own relative coordinate system, where the target tree is at (0, 0) and the neighbor is at (p, q). Now, in the two equations just given, b = 0, h = p, and k = q. In addition, p = x2 - x1 and q = y2 - y1. So now the equations become
Now substitute for y in the second equation, which gives:
If a tree intersects the azimuth angle, it is assessed to see if it also intersects the altitude angle. To determine the altitude, the code again makes the geometry into a two-dimensional problem along the plane of the line projecting from the light ray origin along the azimuth line. In this case the arrangement looks like this:
Δb is the difference between the height of the light ray origin and the base of the neighbor canopy; Δh is the difference between the height of the light ray origin point and the top of the neighbor canopy; dn is the difference between the target tree trunk and the nearest edge of the neighbor canopy; df is the difference between the target tree trunk and the farthest edge of the neighbor canopy. The code needs to solve for θ1 and θ2 to determine if the light ray's altitude angle is between them.
To find dn and df, the code finishes solving the line-and-circle system of equations for the two x solutions - call them x3 and x4. These are plugged into the “y = mx” equation to get the y values - call them y3 and y4. The distance from (0, 0) (in our relative coordinate system, the position of the origin of the light ray) to (x3, y3) is dn; from (0, 0) to (x4, y4) is df. So
After solving for θ1 and θ2, if the light ray's altitude angle is between them, the tree intersects the ray.
In some cases, a tree's canopy extends out above the light ray's origin point. In this case, assessing whether the ray intersects is done a little differently.
In this situation, the canopy always intersects the azimuth angle of the light ray and extends from the altitude angle of the bottom of the canopy (θ) to the zenith. So it only needs to be determined whether the altitude angle of the light ray is above or below the altitude angle of the bottom of the canopy in the plane along the azimuth angle (θ).
The code starts by solving the circle-and-line system of equations as described above to find the two points at which a line of the correct slope intersects the tree canopy. Again, though, we are working with a ray, not a line; one of the points is along the azimuth heading and the other is 180 degrees away. The code determines which solution is along the azimuth heading and throws the other away.
Call the intersection point of the azimuth line and the tree canopy (X3, Y3). Then: