From 1134fee1baccee11cd410ab1db3644bbbc217e29 Mon Sep 17 00:00:00 2001
From: Anton Maminov <anton.linux@gmail.com>
Date: Fri, 29 Mar 2024 12:20:42 +0200
Subject: [PATCH 1/3] add a method to calculate the initial bearing between two
 points

---
 spec/geo/bearing_spec.cr | 18 ++++++++++++++++++
 src/geo/bearing.cr       | 35 +++++++++++++++++++++++++++++++++++
 src/geo/utils.cr         |  8 ++++++++
 3 files changed, 61 insertions(+)
 create mode 100644 spec/geo/bearing_spec.cr
 create mode 100644 src/geo/bearing.cr

diff --git a/spec/geo/bearing_spec.cr b/spec/geo/bearing_spec.cr
new file mode 100644
index 0000000..c7615f4
--- /dev/null
+++ b/spec/geo/bearing_spec.cr
@@ -0,0 +1,18 @@
+require "../spec_helper"
+require "../../src/geo/bearing.cr"
+
+describe Geo::Coord do
+  point_a = Geo::Coord.new(48.8566, 2.3522)   # Paris
+  point_b = Geo::Coord.new(40.7128, -74.0060) # New York
+
+  context "bearing" do
+    # Checked on https://www.movable-type.co.uk/scripts/latlong.html
+    context "calculates initial bearing" do
+      it { point_a.bearing(point_b).should eq(291.7938627483058) }
+    end
+
+    context "calculates final bearing" do
+      it { point_a.bearing(point_b, true).should eq(233.70448129781204) }
+    end
+  end
+end
diff --git a/src/geo/bearing.cr b/src/geo/bearing.cr
new file mode 100644
index 0000000..3798e5c
--- /dev/null
+++ b/src/geo/bearing.cr
@@ -0,0 +1,35 @@
+module Geo
+  struct Coord
+    # The formula used for calculating the initial bearing between two points on the Earth's surfaceis
+    # is derived from the broader concepts of great-circle distance and navigation on a spherical Earth model.
+    # https://en.wikipedia.org/wiki/Great-circle_distance
+    #
+    # https://github.com/Turfjs/turf/blob/master/packages/turf-bearing/index.ts
+    def bearing(to : Geo::Coord, final = false) : Float64
+      if final
+        # Calculate the bearing from the destination point back to the original point
+        reverse_bearing = calculate_bearing(to.lat, to.lng, lat, lng)
+
+        # Adjust by 180 degrees to get the final bearing in the correct direction
+        (reverse_bearing + 180) % 360
+      else
+        calculate_bearing(lat, lng, to.lat, to.lng)
+      end
+    end
+
+    private def calculate_bearing(lat1, lng1, lat2, lng2) : Float64
+      rad_lat1 = Geo::Utils.degrees_to_radians(lat1)
+      rad_lng1 = Geo::Utils.degrees_to_radians(lng1)
+      rad_lat2 = Geo::Utils.degrees_to_radians(lat2)
+      rad_lng2 = Geo::Utils.degrees_to_radians(lng2)
+
+      delta_lng = rad_lng2 - rad_lng1
+
+      a = Math.sin(delta_lng) * Math.cos(rad_lat2)
+      b = Math.cos(rad_lat1) * Math.sin(rad_lat2) -
+          Math.sin(rad_lat1) * Math.cos(rad_lat2) * Math.cos(delta_lng)
+
+      Geo::Utils.radians_to_degrees(Math.atan2(a, b)) % 360
+    end
+  end
+end
diff --git a/src/geo/utils.cr b/src/geo/utils.cr
index 0f27a57..2d7b5be 100644
--- a/src/geo/utils.cr
+++ b/src/geo/utils.cr
@@ -21,5 +21,13 @@ module Geo
       return 0 if val == 0 # colinear
       val > 0 ? 1 : 2      # clockwise or counterclockwise
     end
+
+    def degrees_to_radians(degrees : Number) : Float64
+      degrees * Math::PI / 180.0
+    end
+
+    def radians_to_degrees(radians : Number) : Float64
+      radians * 180.0 / Math::PI
+    end
   end
 end

From e5a79b488c50f5cbb0dd428fd96b3e74bc6490f3 Mon Sep 17 00:00:00 2001
From: Anton Maminov <anton.linux@gmail.com>
Date: Fri, 29 Mar 2024 12:30:52 +0200
Subject: [PATCH 2/3] cosmetic changes

---
 src/geo/polygon.cr | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/src/geo/polygon.cr b/src/geo/polygon.cr
index c711ff4..e952415 100644
--- a/src/geo/polygon.cr
+++ b/src/geo/polygon.cr
@@ -27,7 +27,7 @@ module Geo
 
     # Evaluate area of a polygon using shoelace formula
     def area
-      @coords.each_cons(2).sum { |p| p[0].shoelace(p[1]) }.abs.fdiv(2)
+      @coords.each_cons_pair.sum { |lat, lng| lat.shoelace(lng) }.abs.fdiv(2)
     end
 
     # Order coords in lexicographical order.
@@ -64,17 +64,16 @@ module Geo
       x = coord.lng
       y = coord.lat
 
-      @coords.each do |p|
-        yi = p.lat
-        xi = p.lng
-        yj = last_coord.lat
-        xj = last_coord.lng
+      @coords.each do |iter_coord|
+        yi, xi = iter_coord.ll
+        yj, xj = last_coord.ll
+
         if yi < y && yj >= y ||
            yj < y && yi >= y
           odd_nodes = !odd_nodes if xi + (y - yi) / (yj - yi) * (xj - xi) < x
         end
 
-        last_coord = p
+        last_coord = iter_coord
       end
 
       odd_nodes

From 5fe6e6fd868917779185335f127632f977157a5c Mon Sep 17 00:00:00 2001
From: Anton Maminov <anton.linux@gmail.com>
Date: Fri, 29 Mar 2024 12:33:05 +0200
Subject: [PATCH 3/3] Update polygon.cr

---
 src/geo/polygon.cr | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geo/polygon.cr b/src/geo/polygon.cr
index e952415..39f15ec 100644
--- a/src/geo/polygon.cr
+++ b/src/geo/polygon.cr
@@ -61,8 +61,7 @@ module Geo
     def contains?(coord : Geo::Coord) : Bool
       last_coord = @coords.last
       odd_nodes = false
-      x = coord.lng
-      y = coord.lat
+      y, x = coord.ll
 
       @coords.each do |iter_coord|
         yi, xi = iter_coord.ll
@@ -91,6 +90,7 @@ module Geo
 
     def ==(other : Geo::Polygon) : Bool
       min_size = Math.min(size, other.size)
+
       0.upto(min_size - 1) do |i|
         return false unless self[i] == other[i]
       end