Hurmet function that renders a beam diagram

beam.txt

# Hurmet module for beam diagram creation.
# function diagram() returns a rendered beam diagram.
# function shearMaximums() returns a vector containing the maximum and minimum shear values
# function momentMaximums() returns a vector containing the maximum and minimum bending moment values
# Copyright 2022 Ron Kok
# Released under terms of the MIT License, https://opensource.org/licenses/MIT

function diagram(diagramInput)
   # This function creates a beam diagram.
   # It displays a load diagram, a shear diagram, and a moment diagram.
   # Given a beam divided into n segments, the input to this function is:
   #   (P₁, L₁, w₁ ; P₂, L₂, w₂ ; … ; Pₙ, Lₙ, wₙ), where:
   # P is the point load at the left end of the segment (kips),
   # L is the length of the segment (feet), and
   # w is the line load on the segment (kips/ft).
   # If a point load or reaction exists at the right end, append a [ P, 0, 0 ] to the input matrix.
   
   # Collect information
   lenInput = length(diagramInput) / 3
   numSegments = {
      lenInput     if diagramInput[lenInput, 2] > 0 ;
      lenInput - 1 otherwise
   }
   
   # Load input into vectors
   L = diagramInput[1:numSegments, 2]
   w = diagramInput[1:numSegments, 3]
   P = diagramInput[, 1]
   if lenInput == numSegments
      P = P & 0
   end if
   L_beam = sum(L)
   w_max = max(w)
   w_min = min(w)
   if w_max < 0
      w_max = 0
   end
   if w_min > 0
      w_min = 0
   end
   
   # Calculate shears and append each shear value to a vector.
   # Along the way, keep track of Vmax and Vmin.
   V = [0, P[1]]
   Vmax = 0
   Vmin = 0
   x = 0
   for i in 1:numSegments
      if L[i] < 0
         throw "Error. One segment has a negative length."
      end
      Vstart = V[length(V)]
      if Vstart > Vmax
         Vmax = Vstart
         xVmax = x
      end
      if Vstart < Vmin
         Vmin = Vstart
         xVmin = x
      end if
      Vend = Vstart + w[i] * L[i]
      V = V & Vend
      x = x + L[i]
      if Vend > Vmax
         Vmax = Vend
         xVmax = x
      end
      if Vend < Vmin
         Vmin = Vend
         xVmin = x
      end if
      V = V & (Vend + P[i + 1])
   end
         
   # Calculate local maximums for moment
   Mstart = 0
   Mmax = 0
   Mmin = 0
   x = 0
   for i in 1:numSegments
      Vstart = V[2 * i]
      Vend = V[2 * i + 1]
      if Mstart > Mmax
         Mmax = Mstart
         xMmax = x
      end if
      if Mstart < Mmin
         Mmin = Mstart
         xMmin = x
      end if
      Mend = Mstart + Vstart * L[i] + 0.5 * w[i] * (L[i])²
      if Mend > Mmax
         Mmax = Mend
         xMmax = x + L[i]
      end
      if Mend < Mmin
         Mmin = Mend
         xMmin = x + L[i]
      end
      if sign(Vstart) * sign(Vend) == -1
         xcross = Vstart / -w[i]
         Mmid = Mstart + Vstart * xcross + 0.5 * w[i] * xcross²
        if Mmid > Mmax
           Mmax = Mmid
           xMmax = x + xcross
        end
        if Mmid < Mmin
           Mmin = Mmid
           xMmin = x + xcross
        end
      end
      Mstart = Mend
      x = x + L[i]
   end

   # Work out some geometry
   L_diagram = 300   # width of each diagram, in px
   xScale = L_diagram / L_beam
   wScale = {20 / (w_max - w_min) if (w_max > 0 or w_min < 0); 0 otherwise}
   vScale = 60 / (Vmax - Vmin)
   mScale = 60 / (Mmax - Mmin)
   # y coord of moment diagram is zero
   yV = mScale * Mmax + 40 - vScale * Vmin        # y coord of shear diagram
   yW = yV + vScale * Vmax + 70 - wScale * w_min  # y coord of load diagram
 
   # Now we start to draw things.
   svg = startSvg()
   title "Beam diagram showing loads, shear, amd moment"
   frame 390, 320
   view -65, 315, -30 + {mScale * Mmin if Mmin < 0; 0 otherwise}
   fontsize 12
   
   # Draw baselines for each diagram
   strokewidth 1.5
   line [0, yW; L_diagram, yW]
   text [-60, yW], "loads", "right"
   text [-60, yW - 15], "(kips, ft)", "right"
   strokewidth 1
   line [0, yV; L_diagram, yV]
   text [-60, yV], "shear", "right"
   text [-60, yV - 15], "(kips)", "right"
   # The left end of the moment diagram is set on the origin
   line [0, 0; L_diagram, 0]
   text [-60, 0], "bending", "right"
   text [-60, -15], "(kip·ft)", "right"

   # Draw vertical point loads
   strokewidth 2
   marker "arrow"
   x = 0
   for i in 1:(numSegments + 1)
     if P[i] < -0.005
        line [x * xScale, yW + 30; x * xScale, yW]
        text [x * xScale, yW + 30], string(-P[i], "r3"), "above"
     elseif P[i] > 0.005
        line [x * xScale, yW - 30; x * xScale, yW]
        text [x * xScale, yW - 30], string(P[i], "r3"), "below"
     end
     if i ≤ numSegments
        x = x + L[i]
     end
   end

   # Draw line loads
   strokewidth 1
   fill "#eee"
   marker "none"
   x = 0
   linePath = [0, yW; 0, yW - w[1] * wScale]
   for i in 1:numSegments
      if i > 1
         linePath = vcat(linePath, [x * xScale, yW - w[i] * wScale])
      end
      if L[i] / L_beam > 0.15
         if abs(w[i]) > 0
            pos = {"below" if w[i] > 0; "above" otherwise}
            text [(x + L[i] / 2) * xScale, yW - w[i] * wScale], string(abs(w[i]), "r3"), pos
         end
         pos = {"above" if w[i] > 0; "below" otherwise}
         text [(x + L[i] / 2) * xScale, yW], string(L[i], "r3") & "′", pos
      end if
      x = x + L[i]
      linePath = vcat(linePath, [x * xScale, yW - w[i] * wScale])
   end
   linePath = vcat(linePath, [L_diagram, yW])
   linePath = vcat(linePath, [0, yW])
   path linePath

   # Draw shear diagram
   fill "none"
   x = 0
   shearPath = [0, yV; 0, yV + vScale * V[2]]
   for i in 1:numSegments
      k = 2 * i + 1
      x = x + L[i]
      shearPath = vcat(shearPath, [x * xScale, yV + vScale * V[k]])
      k = k + 1
      shearPath = vcat(shearPath, [x * xScale, yV + vScale * V[k]])
   end
   path shearPath
   if Vmax > 0.001
      text [xVmax * xScale, yV + Vmax * vScale], string(Vmax, "r3"), "above"
   end if
   if Vmin < -0.001
      text [xVmin * xScale, yV + Vmin * vScale], string(-Vmin, "r3"), "below"
   end if

   # Draw moment diagram
   # We'll draw ~75 points
   dx = L_beam / 75
   xStart = 0
   Mstart = 0
   momentPath = [0, 0]
   for i in 1:numSegments
      Vstart = V[2 * i]
      # Vectorize the moment calculation
      x = [0:dx:L[i]...]
      M = (Mstart + Vstart * x + 0.5 * w[i] * (x ∘ x)) * mScale
      x = (x + xStart) * xScale
      momentPath = vcat(momentPath, (x & M))
      xStart = xStart + L[i]
      Mstart = Mstart + Vstart * L[i] + 0.5 * w[i] * (L[i])²
   end
   path momentPath
   if Mmax > 0.001
      text [xMmax * xScale, Mmax * mScale], string(Mmax, "r3"), "above"
   end if
   if Mmin < -0.001
      text [xMmin * xScale, Mmin * mScale], string(-Mmin, "r3"), "below"
   end if
end

function shearMaximums(diagramInput)
   # Collect information
   lenInput = length(diagramInput) / 3
   numSegments = {
      lenInput     if diagramInput[lenInput, 2] > 0 ;
      lenInput - 1 otherwise
   }

   # Load input into vectors
   L = diagramInput[1:numSegments, 2]
   w = diagramInput[1:numSegments, 3]
   P = diagramInput[, 1]
   if lenInput == numSegments
      P = P & 0
   end if

   # Calculate shears
   Vprev = P[1]
   Vmax = 0
   Vmin = 0
   for i in 1:numSegments
      Vstart = Vprev
      Vend = Vstart + w[i] * L[i]
      Vmax = max(Vmax, Vstart, Vend)
      Vmin = min(Vmin, Vstart, Vend)
      Vprev = Vend + P[i + 1]
   end
   return [Vmax, Vmin]
end

function momentMaximums(diagramInput)
   # Collect information
   lenInput = length(diagramInput) / 3
   numSegments = {
      lenInput     if diagramInput[lenInput, 2] > 0 ;
      lenInput - 1 otherwise
   }

   # Load input into vectors
   L = diagramInput[1:numSegments, 2]
   w = diagramInput[1:numSegments, 3]
   P = diagramInput[, 1]
   if lenInput == numSegments
      P = P & 0
   end if

   # Calculate shears
   V = [0, P[1]]
   for i in 1:numSegments
      Vstart = V[length(V)]
      Vend = Vstart + w[i] * L[i]
      V = V & Vend
      V = V & (Vend + P[i + 1])
   end

   # Calculate local maximums for moment
   Mstart = 0
   Mmax = 0
   Mmin = 0
   for i in 1:numSegments
      Vstart = V[2 * i]
      Vend = V[2 * i + 1]
      Mend = Mstart + Vstart * L[i] + 0.5 * w[i] * (L[i])²
      if sign(Vstart) * sign(Vend) == -1
         xcross = Vstart / -w[i]
         Mmid = Mstart + Vstart * xcross + 0.5 * w[i] * xcross²
      else
         Mmid = 0
      end
      Mmax = max(Mmax, Mstart, Mmid, Mend)
      Mmin = min(Mmin, Mstart, Mmid, Mend)
      Mstart = Mend
   end
   return [Mmax, Mmin]
end