From 5310734d4c4c833e025b7140a5eb59a2b9b5155e Mon Sep 17 00:00:00 2001
From: Eatdirt <eatdirt@protonmail.com>
Date: Thu, 20 Feb 2020 21:33:37 +0100
Subject: [PATCH] Add two new methods for geo objects

The method horizon() returns the distance to the horizon, along earth
surface, starting from current location. The method
greatcircle_distance_to(A,B) returns the shortest distance to a great
circle passing by A and B, along earth surface and from current location.
---
 Nasal/geo.nas | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 62 insertions(+)

diff --git a/Nasal/geo.nas b/Nasal/geo.nas
index ac45f970f..18971c6b9 100644
--- a/Nasal/geo.nas
+++ b/Nasal/geo.nas
@@ -45,6 +45,9 @@
 #                                     useful for map distance
 #     .direct_distance_to(<coord>)      ...   distance in m direct, considers altitude,
 #                                             but cuts through Earth surface
+#     .greatcircle_distance_to(<coord>, <coord>)  ... returns distance to a great circle (in m along Earth curvature)
+#                                                     defined by two points  
+#     .horizon()                  ... returns distance to the horizon in m along Earth curvature, ignoring altitudes
 #
 #
 # MANIPULATION METHODS:
@@ -221,6 +224,65 @@ var Coord = {
 		var dz = dest._z - me._z;
 		return math.sqrt(dx * dx + dy * dy + dz * dz);
 	},
+	# arc distance on an earth sphere to the great circle passing by A and B; doesn't consider altitude
+	greatcircle_distance_to: func(destA, destB) {
+		me._pupdate();
+		destA._pupdate();
+		destB._pupdate();
+
+		# AB is not a circle but a point
+		if (destA._lat == destB._lat and destA._lon == destB._lon) {
+		    return me.distance_to(destA);
+		}
+
+		var ca1 = math.cos(destA._lon);
+		var cd1 = math.cos(destA._lat);   
+		var sa1 = math.sin(destA._lon);
+		var sd1 = math.sin(destA._lat);    
+
+		var ca2 = math.cos(destB._lon);
+		var cd2 = math.cos(destB._lat);    
+		var sa2 = math.sin(destB._lon);
+		var sd2 = math.sin(destB._lat);    
+
+		var sa12 = math.sin(destA._lon - destB._lon);
+		
+    		var ca3 = math.cos(me._lon);
+		var cd3 = math.cos(me._lat);    
+		var sa3 = math.sin(me._lon);
+		var sd3 = math.sin(me._lat);    
+
+		# this is sin(greatcircle_dist) * sin(arcAB)
+                var sDsAB = cd3 * sa3 * (ca2 * cd2 * sd1 - ca1 * cd1 * sd2 )
+		    + ca3 * cd3 * ( cd1 * sa1 * sd2 - cd2 * sa2 * sd1 )
+		    - cd1 * cd2 * sd3 * sa12;
+		
+		# direct calculation of sin(arcAB) to not call sin(arcsin(distance_to))    
+		var a = math.sin((destA._lat - destB._lat) * 0.5);
+		var o = math.sin((destA._lon - destB._lon) * 0.5);
+
+		var hs12 = a * a + cd1 * cd2 * o * o;
+		var hc12 = 1.0 - hs12;		
+
+
+		# AB is undertermined; a great circle should be defined with non-colinear vectors
+		if (hs12*hc12 == 0.0) {
+		    die("Great circles are defined with non-colinear vectors");
+		}
+
+		
+		return ERAD * math.abs( math.asin( 0.5 * sDsAB / math.sqrt( hs12 * hc12 ) ) );
+	},
+	# arc distance on an earth sphere to the horizon    
+        horizon: func() {
+	        me._pupdate();
+		if (me._alt < 0.0) {
+		    return 0.0;
+		}
+		else {
+		    return ERAD*math.acos(ERAD/(ERAD+me._alt));
+		}
+	},		
 	is_defined: func {
 		return !(me._cdirty and me._pdirty);
 	},