# HG changeset patch
# User Matti Hamalainen <ccr@tnsp.org>
# Date 1483933540 -7200
# Node ID 3c816fdc302f9e493a528dc0bd635e7c735e4334
# Parent  5edaebbbb7f220dc6a5389134c683773e4d520f2
weather: Implement the "great circle distance" algorithm for calculating
distances between two latitude/longitude coordinate pairs, which should be
more accurate than simply using Euclidean distance.

diff -r 5edaebbbb7f2 -r 3c816fdc302f weather.tcl
--- a/weather.tcl	Mon Jan 09 04:19:10 2017 +0200
+++ b/weather.tcl	Mon Jan 09 05:45:40 2017 +0200
@@ -327,9 +327,9 @@
     }
 
     # Check argument types
-    set d_lat [lindex $nlist 0]
-    set d_lng [lindex $nlist 1]
-    if {![string is double -strict $d_lat] || ![string is double -strict $d_lng]} {
+    set pt1_lat [lindex $nlist 0]
+    set pt1_lng [lindex $nlist 1]
+    if {![string is double -strict $pt1_lat] || ![string is double -strict $pt1_lng]} {
       weather_msg $upublic $unick $uchan $weather_msg_usage_nearest_invalid
       return 0
     }
@@ -338,9 +338,27 @@
     set result {}
     foreach {ukey uvalue} [array get weather_data] {
       if {![string match "w_*" $ukey]} {
-        set delta_lat [expr {$d_lat - [lindex $uvalue 2]}]
-        set delta_lng [expr {$d_lng - [lindex $uvalue 3]}]
-        set dist [expr { sqrt($delta_lat * $delta_lat + $delta_lng * $delta_lng) }]
+        set pt2_lat [lindex $uvalue 2]
+        set pt2_lng [lindex $uvalue 3]
+
+        # Calculate "great circle angle" of two lat/long coordinate pairs
+        # https://en.wikipedia.org/wiki/Great-circle_distance
+        #
+        # The result then can be used to calculate approximate distance
+        # between the points when radius of the sphere "R" is known.
+        set d_lng [expr { $pt2_lng - $pt1_lng }]
+        set c_lat1 [expr { cos($pt1_lat) }]
+        set c_lat2 [expr { cos($pt2_lat) }]
+        set s_lat1 [expr { sin($pt1_lat) }]
+        set s_lat2 [expr { sin($pt2_lat) }]
+
+        set f1 [expr { $c_lat2 * sin($d_lng) }]
+        set f2 [expr { $c_lat2 * cos($d_lng) }]
+        set f3 [expr { $c_lat1 * $s_lat2 - $s_lat1 * $f2 }]
+        set c_ang [expr { atan2(sqrt($f1 * $f1 + $f3 * $f3), $s_lat1 * $s_lat2 + $c_lat1 * $f2) }]
+
+        # Calculate distance based on Earth radius approximation in km * 1000m
+        set dist [expr { 6371.439 * 1000.0 * $c_ang }]
         lappend result [list $ukey $dist]
       }
     }
@@ -356,7 +374,7 @@
 
     # Print out the result
     set res [join $uresult " ; "]
-    weather_msg $upublic $unick $uchan $weather_msg_nearest_stations [list $d_lat $d_lng $res]
+    weather_msg $upublic $unick $uchan $weather_msg_nearest_stations [list $pt1_lat $pt1_lng $res]
     return 0
   } elseif {$rarg == "vakio" || $rarg == "default" || $rarg == "vakiot" || $rarg == "defaults"} {
     # List or set the default weather station name patterns for this user