Calculate Growing Degree Days (GDD) in ruby

gdd.rb

# Calculates Growing Degree Days (GDD) using generally-accepted methods
# defined by the Division of Agriculture and Natural Resources, University
# of California [1].
#
# The single-sine method seems to be the defacto way of determining a GDD [2].
#
# The triangulation methods are originally from [3] and are used to correct
# a transcription error in [1].
#
#
# References:
#
# [1] Degree-days: the Calculation And Use of Heat Units In Pest Management. Berkeley, Calif.:
#     Division of Agriculture and Natural Resources, University of California, 1983.
#
# [2] http://ipm.ucanr.edu/WEATHER/ddeval.html
#
# [3] Heat accumulation for timing Lygus control measures in a safflower-cotton complex:
#     V. Sevacherian, V. M. Stern, A. J. Mueller, 1977.
#
class GrowingDegreeDay

  def self.double_triangle(tmin1, tmax, tmin2, tfloor, tceil)
    half_triangle(tmin1, tmax, tfloor, tceil) +
    half_triangle(tmin2, tmax, tfloor, tceil)
  end

  def self.single_triangle(tmin, tmax, tfloor, tceil)
    half_triangle(tmin, tmax, tfloor, tceil) * 2
  end

  def self.double_sine(tmin1, tmax, tmin2, tfloor, tceil)
    half_sine(tmin1, tmax, tfloor, tceil) +
    half_sine(tmin2, tmax, tfloor, tceil)
  end

  def self.single_sine(tmin, tmax, tfloor, tceil)
    half_sine(tmin, tmax, tfloor, tceil) * 2
  end

  def self.half_triangle(tmin, tmax, tfloor, tceil)
    case bounding(tmin, tmax, tfloor, tceil)

    when BOUND_ABOVE_BOTH
      (tceil - tfloor) / 2.0

    when BOUND_BELOW_BOTH
      0

    when BOUND_BETWEEN_BOTH
      (6 * (tmax + tmin - (2 * tfloor))) / 24.0

    when BOUND_INTERCEPT_LOWER
      tdiff = (tmax - tmin).to_f

      # The formula as transcribed in the original 1983 leaflet [1] is incorrect.
      # The correct version is used here. The updated version is also available
      # at http://ipm.ucanr.edu/WEATHER/ddfig44.html.

      6 * (tmax - tfloor) ** 2 / tdiff / 24.0

    when BOUND_INTERCEPT_UPPER
      tdiff = (tmax - tmin).to_f

      x = 6 * (tmax + tmin - 2 * tfloor) / 24.0
      y = 6 * (tmax - tceil) ** 2 / tdiff

      x - y / 24.0


    when BOUND_INTERCEPT_BOTH
      tdiff = (tmax - tmin).to_f

      x = 6 * (tmax - tfloor) ** 2 / tdiff
      y = 6 * (tmax - tceil) ** 2 / tdiff

      (x - y) / 24.0

    end
  end

  # The period of the sine curve.
  TWO_PI  = Math::PI * 2
  HALF_PI = Math::PI / 2

  def self.half_sine(tmin, tmax, tfloor, tceil)
    case bounding(tmin, tmax, tfloor, tceil)

    when BOUND_ABOVE_BOTH
      (tceil - tfloor) / 2.0

    when BOUND_BELOW_BOTH
      0

    when BOUND_BETWEEN_BOTH
      0.5 * ((tmax + tmin) / 2.0 - tfloor)

    when BOUND_INTERCEPT_LOWER
      a = (tmax - tmin) / 2.0
      x = (tmax + tmin) / 2.0

      o1 = Math::asin( (tfloor - x) / a )

      1/TWO_PI * ( (x - tfloor) * (HALF_PI - o1) + a * Math::cos(o1) )

    when BOUND_INTERCEPT_UPPER
      a = (tmax - tmin) / 2.0
      x = (tmax + tmin) / 2.0

      o2 = Math::asin( (tceil - x) / a )

      1/TWO_PI * ( (x-tfloor) * (o2 + HALF_PI) + (tceil - tfloor) * (HALF_PI - o2) - (a * Math::cos(o2) ) )

    when BOUND_INTERCEPT_BOTH
      a = (tmax - tmin) / 2.0
      x = (tmax + tmin) / 2.0

      o1 = Math::asin( (tfloor - x) / a )
      o2 = Math::asin( (tceil - x) / a )

      1/TWO_PI * ( (x-tfloor) * (o2 - o1) + ( a * (Math::cos(o1) - Math::cos(o2)) ) + (tceil - tfloor) * (HALF_PI - o2) )
    end
  end

  # different bounds, referred as cases 1-6 in [1].

  # Above both thresholds.
  BOUND_ABOVE_BOTH = :aboveboth

  # Below both thresholds.
  BOUND_BELOW_BOTH = :belowboth

  # Between both thresholds.
  BOUND_BETWEEN_BOTH = :betweenboth

  # Intercepted by the lower threshold.
  BOUND_INTERCEPT_LOWER = :interlower

  # Intercepted by the upper threshold.
  BOUND_INTERCEPT_UPPER = :interupper

  # Intercepted by both thresholds.
  BOUND_INTERCEPT_BOTH = :interboth

  def self.bounding(tmin, tmax, tfloor, tceil)
    thresholds = tfloor..tceil
    temps = tmin..tmax

    case
    when tmin > tceil then BOUND_ABOVE_BOTH
    when tmax < tfloor then BOUND_BELOW_BOTH
    when thresholds.cover?(temps) then BOUND_BETWEEN_BOTH
    when thresholds.include?(tmax) && !thresholds.include?(tmin) then BOUND_INTERCEPT_LOWER
    when thresholds.include?(tmin) && !thresholds.include?(tmax) then BOUND_INTERCEPT_UPPER
    when temps.cover?(thresholds) then BOUND_INTERCEPT_BOTH
    else raise "Bounds calculation error"
    end
  end

end

test_gdd.rb

require "gdd"

# GDD test cases from https://beaumont.tamu.edu/eLibrary/Publications/Ted_Wilson/LTW31.pdf page 9.
# https://catalog.hathitrust.org/Record/008707238. Some values have been corrected using the
# GDD calculator at http://ipm.ucanr.edu/WEATHER/ddretrieve.html. They are marked with an asterisk.
#
# Transcription notes:
#
#  The reference in the text to `sin -1` should be read as the `arcsin` function. See
#  https://en.wikipedia.org/w/index.php?title=Inverse_trigonometric_functions&oldid=953132433#Notation
#
# The temperature cycle can be:
#  1) completely above both thresholds,
#  2) completely below both thresholds,
#  3) entirely between both thresholds,
#  4) intercepted by the lower threshold,
#  5) intercepted by the upper threshold, or
#  6) intercepted by both thresholds.
#
# Development Thresholds: 55F min, 90F max
#
# Cycle Min  Max  Min    Triangulation        Sine
# Type  temp temp temp   Single  Double   Single  Double
#
#   1    96   110  91     35      35       35      35
#   1    91   105  96     35      35       35      35
#
#   2    45   54   38      0       0        0       0
#   2    38   50   45      0       0        0       0
#
#   3    60   80   75     15.00   18.75    15.00   18.75
#   3    75   88   60     26.50   22.75    26.50   22.75
#
#   4    50   82   45     11.39   10.62    11.85   11.31
#   4    45   70   50      4.50    5.06     5.31   5.70
#
#   5    60   100  75     23.75   27.13*   22.82   26.26
#   5    75   95   60     29.38   25.76    28.91   25.30
#
#   6    50   101  48     19.56   19.19    18.95   18.69
#   6    48   95   50     16.76   17.13    16.96*  17.23
#

require 'minitest/autorun'

class TestGdd < Minitest::Test

  # Thresholds
  TLO = 55
  THI = 90

  def test_above_both
    assert_in_delta 35.0, GrowingDegreeDay.single_triangle(96,110,TLO,THI)
    assert_in_delta 35.0, GrowingDegreeDay.single_triangle(91,105,TLO,THI)

    assert_in_delta 35.0, GrowingDegreeDay.single_sine(96,110,TLO,THI)
    assert_in_delta 35.0, GrowingDegreeDay.single_sine(91,105,TLO,THI)

    assert_in_delta 35.0, GrowingDegreeDay.double_triangle(96,110,91,TLO,THI)
    assert_in_delta 35.0, GrowingDegreeDay.double_triangle(91,105,96,TLO,THI)

    assert_in_delta 35.0, GrowingDegreeDay.double_sine(96,110,91,TLO,THI)
    assert_in_delta 35.0, GrowingDegreeDay.double_sine(91,105,96,TLO,THI)
  end

  def test_below_both
    assert_equal 0, GrowingDegreeDay.single_triangle(45,54,TLO,THI)
    assert_equal 0, GrowingDegreeDay.single_triangle(38,50,TLO,THI)

    assert_equal 0, GrowingDegreeDay.single_sine(45,54,TLO,THI)
    assert_equal 0, GrowingDegreeDay.single_sine(38,50,TLO,THI)

    assert_equal 0, GrowingDegreeDay.double_triangle(45,54,38,TLO,THI)
    assert_equal 0, GrowingDegreeDay.double_triangle(38,50,45,TLO,THI)

    assert_equal 0, GrowingDegreeDay.double_sine(45,54,38,TLO,THI)
    assert_equal 0, GrowingDegreeDay.double_sine(38,50,45,TLO,THI)
  end

  def test_between_both
    assert_in_delta 15.0, GrowingDegreeDay.single_triangle(60,80,TLO,THI)
    assert_in_delta 26.5, GrowingDegreeDay.single_triangle(75,88,TLO,THI)

    assert_in_delta 18.75, GrowingDegreeDay.double_triangle(60,80,75,TLO,THI)
    assert_in_delta 22.75, GrowingDegreeDay.double_triangle(75,88,60,TLO,THI)

    assert_in_delta 15.0, GrowingDegreeDay.single_sine(60,80,TLO,THI)
    assert_in_delta 26.5, GrowingDegreeDay.single_sine(75,88,TLO,THI)

    assert_in_delta 18.75, GrowingDegreeDay.double_sine(60,80,75,TLO,THI)
    assert_in_delta 22.75, GrowingDegreeDay.double_sine(75,88,60,TLO,THI)
  end

  def test_intercept_lower
    assert_in_delta 11.39, GrowingDegreeDay.single_triangle(50,82,TLO,THI), 0.01
    assert_in_delta  4.5,  GrowingDegreeDay.single_triangle(45,70,TLO,THI), 0.01

    assert_in_delta 10.62, GrowingDegreeDay.double_triangle(50,82,45,TLO,THI), 0.01
    assert_in_delta  5.06, GrowingDegreeDay.double_triangle(45,70,50,TLO,THI), 0.01

    assert_in_delta 11.85, GrowingDegreeDay.single_sine(50,82,TLO,THI), 0.01
    assert_in_delta  5.31, GrowingDegreeDay.single_sine(45,70,TLO,THI), 0.01

    assert_in_delta 11.31, GrowingDegreeDay.double_sine(50,82,45,TLO,THI), 0.01
    assert_in_delta  5.7,  GrowingDegreeDay.double_sine(45,70,50,TLO,THI), 0.01
  end

  def test_intercept_upper
    assert_in_delta 23.75, GrowingDegreeDay.single_triangle(60,100,TLO,THI), 0.01
    assert_in_delta 29.38, GrowingDegreeDay.single_triangle(75,95,TLO,THI), 0.01

    assert_in_delta 27.13, GrowingDegreeDay.double_triangle(60,100,75,TLO,THI), 0.01
    assert_in_delta 25.76, GrowingDegreeDay.double_triangle(75,95,60,TLO,THI), 0.01

    assert_in_delta 22.82, GrowingDegreeDay.single_sine(60,100,TLO,THI), 0.01
    assert_in_delta 28.91, GrowingDegreeDay.single_sine(75,95,TLO,THI), 0.01

    assert_in_delta 26.26, GrowingDegreeDay.double_sine(60,100,75,TLO,THI), 0.01
    assert_in_delta 25.30, GrowingDegreeDay.double_sine(75,95,60,TLO,THI), 0.01
  end

  def test_intercept_both
    assert_in_delta 19.56, GrowingDegreeDay.single_triangle(50,101,TLO,THI), 0.01
    assert_in_delta 16.76, GrowingDegreeDay.single_triangle(48,95,TLO,THI), 0.01

    assert_in_delta 19.19, GrowingDegreeDay.double_triangle(50,101,48,TLO,THI), 0.01
    assert_in_delta 17.13, GrowingDegreeDay.double_triangle(48,95,50,TLO,THI), 0.01

    assert_in_delta 18.95, GrowingDegreeDay.single_sine(50,101,TLO,THI), 0.01
    assert_in_delta 16.96, GrowingDegreeDay.single_sine(48,95,TLO,THI), 0.01

    assert_in_delta 18.69, GrowingDegreeDay.double_sine(50,101,48,TLO,THI), 0.01
    assert_in_delta 17.23, GrowingDegreeDay.double_sine(48,95,50,TLO,THI), 0.01
  end

end

$ ruby -wI. test_gdd.rb -p 
Run options: -p --seed 10332

# Running:

......

Fabulous run in 0.001247s, 4811.5479 runs/s, 38492.3833 assertions/s.

6 runs, 48 assertions, 0 failures, 0 errors, 0 skips