Add interpolation method (#52)

Interpolation is very efficient and useful to compute intermediate
values for instants or quantities that are equidistant.

This adds an interpolation method for three of five terms, given an
interpolation factor.

lib/astronoby.rb

@@ -21,6 +21,7 @@
 require "astronoby/time/greenwich_sidereal_time"
 require "astronoby/time/local_sidereal_time"
 require "astronoby/util/astrodynamics"
+require "astronoby/util/maths"
 require "astronoby/util/time"
 require "astronoby/util/trigonometry"
 require "astronoby/true_obliquity"

lib/astronoby/util/maths.rb

@@ -0,0 +1,72 @@
+# frozen_string_literal: true
+
+module Astronoby
+  module Util
+    module Maths
+      class << self
+        # Source:
+        #  Title: Astronomical Algorithms
+        #  Author: Jean Meeus
+        #  Edition: 2nd edition
+        #  Chapter: 3 - Interpolation
+
+        # @param values [Array<Numeric>] First term
+        # @param factor [Numeric] Interpolation factor
+        # @return [Float] Interpolated value
+        def interpolate(values, factor)
+          unless factor.between?(0, 1)
+            raise IncompatibleArgumentsError,
+              "Interpolation factor must be between 0 and 1"
+          end
+
+          if values.length == 3
+            return interpolate_3_terms(values, factor)
+          elsif values.length == 5
+            return interpolate_5_terms(values, factor)
+          end
+
+          raise IncompatibleArgumentsError,
+            "Only 3 or 5 terms are supported for interpolation"
+        end
+
+        private
+
+        # @return [Float] Interpolated value
+        def interpolate_3_terms(terms, factor)
+          y1, y2, y3 = terms
+
+          a = y2 - y1
+          b = y3 - y2
+          c = b - a
+
+          y2 + (factor / 2.0) * (a + b + factor * c)
+        end
+
+        # @return [Float] Interpolated value
+        def interpolate_5_terms(terms, factor)
+          y1, y2, y3, y4, y5 = terms
+
+          a = y2 - y1
+          b = y3 - y2
+          c = y4 - y3
+          d = y5 - y4
+
+          e = b - a
+          f = c - b
+          g = d - c
+
+          h = f - e
+          j = g - f
+
+          k = j - h
+
+          y3 +
+            factor * ((b + c) / 2.0 - (h + j) / 12.0) +
+            factor**2 * (f / 2.0 - k / 24.0) +
+            factor**3 * (h + j) / 12.0 +
+            factor**4 * k / 24.0
+        end
+      end
+    end
+  end
+end

spec/astronoby/util/maths_spec.rb

@@ -0,0 +1,64 @@
+# frozen_string_literal: true
+
+RSpec.describe Astronoby::Util::Maths do
+  describe "::interpolate" do
+    # Source:
+    #  Title: Astronomical Algorithms
+    #  Author: Jean Meeus
+    #  Edition: 2nd edition
+    #  Chapter: 3 - Interpolation, p.25
+    it "calculates the interpolation for 3 terms" do
+      interpolation = described_class.interpolate(
+        [0.884226, 0.877366, 0.870531],
+        0.18125
+      )
+
+      expect(interpolation.round(6)).to eq 0.876125
+    end
+
+    # Source:
+    #  Title: Astronomical Algorithms
+    #  Author: Jean Meeus
+    #  Edition: 2nd edition
+    #  Chapter: 3 - Interpolation, p.25
+    it "calculates the interpolation for 3 terms" do
+      interpolation = described_class.interpolate(
+        [1.3814294, 1.3812213, 1.3812453],
+        0.3966
+      )
+
+      expect(interpolation.round(6)).to eq 1.381203
+    end
+
+    # Source:
+    #  Title: Astronomical Algorithms
+    #  Author: Jean Meeus
+    #  Edition: 2nd edition
+    #  Chapter: 3 - Interpolation, p.25
+    it "calculates the interpolation for 5 terms" do
+      interpolation = described_class.interpolate(
+        [
+          Astronoby::Angle.from_hms(0, 54, 36.125).hours * 3600,
+          Astronoby::Angle.from_hms(0, 54, 24.606).hours * 3600,
+          Astronoby::Angle.from_hms(0, 54, 15.486).hours * 3600,
+          Astronoby::Angle.from_hms(0, 54, 8.694).hours * 3600,
+          Astronoby::Angle.from_hms(0, 54, 4.133).hours * 3600
+        ],
+        0.27777778
+      )
+
+      expect(interpolation.round(3))
+        .to eq(Astronoby::Angle.from_hms(0, 54, 13.369).hours * 3600)
+    end
+
+    context "when the interpolation factor is out of range" do
+      it "raises an error" do
+        expect { described_class.interpolate([1, 2, 3], 4) }
+          .to raise_error(
+            Astronoby::IncompatibleArgumentsError,
+            "Interpolation factor must be between 0 and 1"
+          )
+      end
+    end
+  end
+end