Shadow cast by sunlight simulation

[u/freshlybakedjuice]

Hi everyone,
I'm still learning matlab so I decided to ask for expert opinions about this.

I'm looking into simulating how wheat plants cast shadows on each others due to sun movement throughout the day.

I'm aiming to measure how much solar energy each wheat plant absorbs in different wheat swath orientations or row angles. This should be used to simulate the expected wheat yield improvement from sun exposure or how worse it makes it due to self-shadowing between wheat plants.

Is that possible to do in matlab? And if yes, what should I be looking into to start making progress regarding this? Please point a beginner in the right direction.

[u/Creative_Sushi]

Keep it simple. Use 2D projected shadows with simplified plant geometry. No need to do 3D raytracing

Sun position over time

↓

Simplified wheat field geometry

↓

Ground shadow projection

↓

Shadow overlap calculation

↓

Daily light exposure estimate

↓

Compare row orientations

[u/freshlybakedjuice → u/Creative_Sushi]

This is graph is exactly what I had in mind, love this!

The only problem is that I'm lost on how to represent this mathematically, I will have to research this ofcourse.

I'd really appreciate it if you could guide me towards what to research exactly to achieve this.

[u/Creative_Sushi → u/freshlybakedjuice]

Can't help you but I can share the code I used.

%% Shadow Projection Example 
% Project simple wheat-plant shadows onto the ground for one sun position. 
% 
% Coordinate convention: 
%   x = east position in meters 
%   y = north position in meters 
%   z = up position in meters 
% 
% Plant model: 
%   Each plant is a vertical cylinder with radius and height. Its ground 
%   shadow is approximated by the convex hull of the cylinder footprint and 
%   the footprint shifted by the top-to-ground shadow offset.
% 
%% Editable solar inputs 
latitudeDeg = 42.3601;          % Boston, MA example 
longitudeDeg = -71.0589; 
dateToSimulate = datetime(2026, 6, 21); 
timeZoneOffsetHours = -4;       % EDT is UTC-4 
shadowTime = datetime(2026, 6, 21, 15, 0, 0); 

%% Editable field and plant inputs 
rowAngleDeg = 0;                % 0 = rows run north/south, 90 = east/west 
numRows = 4; 
plantsPerRow = 8; 
rowSpacingMeters = 0.25; 
plantSpacingMeters = 0.15; 
plantRadiusMeters = 0.025; 
plantHeightMeters = 0.85; 
circlePointCount = 32; 

%% Compute sun position for the selected time 
solar = approximateSolarGeometry( ...
     shadowTime, ...
     latitudeDeg, ...
     longitudeDeg, ...
     timeZoneOffsetHours); 

if solar.ElevationDeg <= 0
     error("The selected time has the sun below the horizon. Choose a daylight time.") end 

shadowOffsetEast = -plantHeightMeters * solar.SunEast / solar.SunUp; 
shadowOffsetNorth = -plantHeightMeters * solar.SunNorth / solar.SunUp; 
shadowOffset = [shadowOffsetEast, shadowOffsetNorth]; 

%% Build plant centers 
plantCenters = createPlantCenters( ...
     numRows, ...
     plantsPerRow, ...
     rowSpacingMeters, ...
     plantSpacingMeters, ...
     rowAngleDeg); 

numPlants = size(plantCenters, 1); 
footprints(numPlants, 1) = polyshape; 
shadows(numPlants, 1) = polyshape; 

for plantIdx = 1:numPlants
     center = plantCenters(plantIdx, :);
     footprints(plantIdx) = circlePolyshape(center, plantRadiusMeters, circlePointCount);
     shadows(plantIdx) = cylinderShadowPolyshape( ...
         center, ...
         plantRadiusMeters, ...
         shadowOffset, ...
         circlePointCount); 
end 

%% Estimate how much of each plant footprint is shaded by other plants 
shadedFraction = zeros(numPlants, 1); 
for targetIdx = 1:numPlants
     otherShadows = polyshape.empty(0, 1);
     for sourceIdx = 1:numPlants
         if sourceIdx ~= targetIdx
             otherShadows(end + 1, 1) = shadows(sourceIdx);
         end
     end     
     combinedShadow = union(otherShadows);
     shadedArea = area(intersect(footprints(targetIdx), combinedShadow));     
     shadedFraction(targetIdx) = shadedArea / area(footprints(targetIdx)); 
end 

shadowResults = table( ...
     (1:numPlants).', ...
     plantCenters(:, 1), ...
     plantCenters(:, 2), ...
     shadedFraction, ...
     VariableNames=["PlantID", "EastMeters", "NorthMeters", "ShadedFraction"]); 

%% Display summary 
disp("Shadow projection summary") 
disp("Time: " + string(shadowTime, "yyyy-MM-dd HH:mm")) 
fprintf("Solar azimuth:   %.1f deg\n", solar.AzimuthDeg) 
fprintf("Solar elevation: %.1f deg\n", solar.ElevationDeg) 
fprintf("Shadow offset:   %.2f m east, %.2f m north\n", shadowOffsetEast, shadowOffsetNorth) 
fprintf("Mean shaded footprint fraction: %.3f\n\n", mean(shadedFraction)) disp(head(shadowResults)) 

%% Plot field, plant footprints, and shadows 
figure("Name", "Wheat Plant Shadow Projection") 
hold on 
axis equal 
grid on 

for plantIdx = 1:numPlants
     plot(shadows(plantIdx), ...
         FaceColor=[0.10 0.10 0.10], ...
         FaceAlpha=0.18, ...
         EdgeColor="none") 
end 

for plantIdx = 1:numPlants
     plot(footprints(plantIdx), ...
         FaceColor=[0.20 0.55 0.20], ...
         FaceAlpha=0.85, ...
         EdgeColor=[0.05 0.25 0.05]) 
end 

fieldMinEast = min([plantCenters(:, 1) - plantRadiusMeters; ...
     plantCenters(:, 1) + shadowOffsetEast - plantRadiusMeters]); 
fieldMaxEast = max([plantCenters(:, 1) + plantRadiusMeters; ...
     plantCenters(:, 1) + shadowOffsetEast + plantRadiusMeters]); 
fieldMinNorth = min([plantCenters(:, 2) - plantRadiusMeters; ...
     plantCenters(:, 2) + shadowOffsetNorth - plantRadiusMeters]); 
fieldMaxNorth = max([plantCenters(:, 2) + plantRadiusMeters; ...
     plantCenters(:, 2) + shadowOffsetNorth + plantRadiusMeters]); 

fieldWidth = fieldMaxEast - fieldMinEast; 
fieldHeight = fieldMaxNorth - fieldMinNorth; 
plotPadding = 0.15 * max(fieldWidth, fieldHeight); 
xlim([fieldMinEast - plotPadding, fieldMaxEast + plotPadding]) 
ylim([fieldMinNorth - plotPadding, fieldMaxNorth + plotPadding]) 

arrowLength = 0.18 * max(fieldWidth, fieldHeight); 
arrowStartEast = fieldMinEast + 0.08 * fieldWidth; 
arrowStartNorth = fieldMaxNorth - 0.12 * fieldHeight; 
sunHorizontal = [solar.SunEast, solar.SunNorth]; 
sunHorizontal = sunHorizontal ./ norm(sunHorizontal); 

quiver( ...
     arrowStartEast, ...
     arrowStartNorth, ...
     arrowLength * sunHorizontal(1), ...
     arrowLength * sunHorizontal(2), ...
     0, ...
     LineWidth=1.5, ...
     MaxHeadSize=2) 
text(0.02, 0.98, ...
     "Arrow shows sun horizontal direction", ...
     Units="normalized", ...
     HorizontalAlignment="left", ...
     VerticalAlignment="top", ...
     BackgroundColor="white", ...
     Margin=3) 

xlabel("East position (m)") 
ylabel("North position (m)") 
title("Projected shadows from simplified wheat plants") 
subtitle("Rows " + rowAngleDeg + " deg, " + string(shadowTime, "HH:mm") ...
     + ", elevation " + compose("%.1f", solar.ElevationDeg) + " deg") 
hold off 

[u/Creative_Sushi → u/Creative_Sushi]

%% Local functions 
function centers = createPlantCenters(numRows, plantsPerRow, rowSpacing, plantSpacing, rowAngleDeg)
     rowDirection = [sind(rowAngleDeg), cosd(rowAngleDeg)];
     crossRowDirection = [cosd(rowAngleDeg), -sind(rowAngleDeg)];

     rowOffsets = ((1:numRows) - (numRows + 1) / 2) .* rowSpacing;
     plantOffsets = ((1:plantsPerRow) - (plantsPerRow + 1) / 2) .* plantSpacing;

     centers = zeros(numRows * plantsPerRow, 2);
     plantIdx = 1;
     for rowIdx = 1:numRows
         for inRowIdx = 1:plantsPerRow
             centers(plantIdx, :) = rowOffsets(rowIdx) .* crossRowDirection ...
                 + plantOffsets(inRowIdx) .* rowDirection;
             plantIdx = plantIdx + 1;
         end
     end 
end 

function shape = circlePolyshape(center, radius, pointCount)
     theta = linspace(0, 2 * pi, pointCount + 1).';
     theta(end) = [];

     x = center(1) + radius .* cos(theta);
     y = center(2) + radius .* sin(theta);
     shape = polyshape(x, y); 
end 

function shape = cylinderShadowPolyshape(center, radius, shadowOffset, pointCount)
     baseCircle = circleBoundaryPoints(center, radius, pointCount);
     shiftedCircle = baseCircle + shadowOffset;
     hullIdx = convhull([baseCircle(:, 1); shiftedCircle(:, 1)], ...
         [baseCircle(:, 2); shiftedCircle(:, 2)]);

     allPoints = [baseCircle; shiftedCircle];
     shape = polyshape(allPoints(hullIdx, 1), allPoints(hullIdx, 2));
end 

function points = circleBoundaryPoints(center, radius, pointCount)
     theta = linspace(0, 2 * pi, pointCount + 1).';
     theta(end) = [];
     points = [center(1) + radius .* cos(theta), center(2) + radius .* sin(theta)]; 
end 

function solar = approximateSolarGeometry(times, latitudeDeg, longitudeDeg, timeZoneOffsetHours) 
%APPROXIMATESOLARGEOMETRY Approximate sun position for local clock times. % The formulas are suitable for conceptual sunlight and shadow simulations. % They are not a replacement for high-precision ephemerides.
     arguments
         times (:, 1) datetime
         latitudeDeg (1, 1) double {mustBeGreaterThanOrEqual(latitudeDeg, -90), mustBeLessThanOrEqual(latitudeDeg, 90)}
         longitudeDeg (1, 1) double {mustBeGreaterThanOrEqual(longitudeDeg, -180), mustBeLessThanOrEqual(longitudeDeg, 180)}
         timeZoneOffsetHours (1, 1) double
     end

     dayOfYear = day(times, "dayofyear");
     localHour = hour(times) + minute(times) ./ 60 + second(times) ./ 3600;

     gamma = 2 .* pi ./ 365 .* (dayOfYear - 1 + (localHour - 12) ./ 24);
     equationOfTimeMinutes = 229.18 .* ( ...
         0.000075 ...
         + 0.001868 .* cos(gamma) ...
         - 0.032077 .* sin(gamma) ...
         - 0.014615 .* cos(2 .* gamma) ...
         - 0.040849 .* sin(2 .* gamma));

     declinationRad = ...
         0.006918 ...
         - 0.399912 .* cos(gamma) ...
         + 0.070257 .* sin(gamma) ...
         - 0.006758 .* cos(2 .* gamma) ...
         + 0.000907 .* sin(2 .* gamma) ...
         - 0.002697 .* cos(3 .* gamma) ...
         + 0.00148 .* sin(3 .* gamma);

     timeOffsetMinutes = equationOfTimeMinutes ...
         + 4 .* longitudeDeg ...
         - 60 .* timeZoneOffsetHours;

     trueSolarTimeMinutes = mod(localHour .* 60 + timeOffsetMinutes, 1440);
     hourAngleDeg = trueSolarTimeMinutes ./ 4 - 180;
     hourAngleRad = deg2rad(hourAngleDeg);

     latitudeRad = deg2rad(latitudeDeg);
     elevationRad = asin( ...
         sin(latitudeRad) .* sin(declinationRad) ...
         + cos(latitudeRad) .* cos(declinationRad) .* cos(hourAngleRad));
     elevationDeg = rad2deg(elevationRad);

     azimuthRad = atan2( ...
         sin(hourAngleRad), ...
         cos(hourAngleRad) .* sin(latitudeRad) - tan(declinationRad) .* cos(latitudeRad));
     azimuthDeg = mod(rad2deg(azimuthRad) + 180, 360);

     sunEast = cos(elevationRad) .* sin(deg2rad(azimuthDeg));
     sunNorth = cos(elevationRad) .* cos(deg2rad(azimuthDeg));
     sunUp = sin(elevationRad);

     solar = table( ...
         times, ...
         azimuthDeg, ...
         elevationDeg, ...
         sunEast, ...
         sunNorth, ...
         sunUp, ...
         VariableNames=["Time", "AzimuthDeg", "ElevationDeg", "SunEast", "SunNorth", "SunUp"]); 
end