From 31ee995825e3a5bc640d747be57a9be109dd060a Mon Sep 17 00:00:00 2001
From: Matt Birkholz <matt@birkholz.chandler.az.us>
Date: Sun, 3 Nov 2013 17:20:47 -0700
Subject: [PATCH] planetarium: Add geodesic-distance, angular-separation,...

...and an angular-separation test procedure.
---
 src/planetarium/geometry.scm | 79 +++++++++++++++++++++++++++++++++++-
 src/planetarium/mit.pkg      |  4 +-
 src/planetarium/time.scm     |  5 ---
 3 files changed, 80 insertions(+), 8 deletions(-)

diff --git a/src/planetarium/geometry.scm b/src/planetarium/geometry.scm
index 45a9f90dc..9026b05c7 100644
--- a/src/planetarium/geometry.scm
+++ b/src/planetarium/geometry.scm
@@ -107,6 +107,78 @@ USA.
 	    (flo:+ r+frac 360.)
 	    r+frac)))))
 
+(define (geodesic-distance p1 p2)
+  ;; "[Given] the geographic coordinates of two points on the surface
+  ;; of the Earth [...] the shortest distance S between these points,
+  ;; measured along the Earth's surface. [...]  We will suppose that
+  ;; these points are at sea level [and] no great accuracy is needed."
+  ;; -Meeus 2009, specifically the second English edition (1998) "with
+  ;; corrections as of August 10, 2009" of _Astronomical_Algorithms_
+  ;; by Jean Meeus.
+
+  (let ((d (angular-separation p1 p2)))
+    (flo:* d (flo:/ (flo:* 6371000.0 pi)
+		    180.0))))
+
+(define (angular-separation p1 p2)
+  ;; "Mr. Thierry Pauwels, of the Royal Observatory of Belgium,
+  ;; communicated the following method. [...]  Mathematically
+  ;; speaking, this method is completely identical to formula 17.1,
+  ;; but a computer will yield more accurate results from an
+  ;; arctangent than from an arccosine." [p.115]
+  ;; -Meeus 2009, specifically the second English edition (1998) "with
+  ;; corrections as of August 10, 2009" of _Astronomical_Algorithms_
+  ;; by Jean Meeus.
+
+  (let ((lat1 (degrees->radians (latitude p1)))
+	(lng1 (degrees->radians (longitude p1)))
+	(lat2 (degrees->radians (latitude p2)))
+	(lng2 (degrees->radians (longitude p2))))
+    (let ((x (flo:- (flo:* (flo:cos lat1) (flo:sin lat2))
+		    (flo:* (flo:sin lat1)
+			   (flo:* (flo:cos lat2)
+				  (flo:cos (flo:- lng2 lng1))))))
+	  (y (flo:* (flo:cos lat2) (flo:sin (flo:- lng2 lng1))))
+	  (z (flo:+ (flo:* (flo:sin lat1) (flo:sin lat2))
+		    (flo:* (flo:cos lat1)
+			   (flo:* (flo:cos lat2)
+				  (flo:cos (flo:- lng2 lng1)))))))
+      (let ((d (flo:atan (flo:/ (flo:sqrt (flo:+ (flo:* x x) (flo:* y y)))
+				z))))
+	(radians->degrees
+	 (if (flo:negative? z)
+	     (flo:+ pi d)
+	     d))))))
+
+(define (test-angular-separation)
+
+  (define (test lat1 lng1 lat2 lng2 degrees tolerance)
+    (let ((d (angular-separation
+	      (make-latitude/longitude lat1 lng1)
+	      (make-latitude/longitude lat2 lng2))))
+      (if (not (flo:~= degrees d tolerance))
+	  (warn "Wrong angular separation:" d degrees))))
+
+  ;; Arcturus to Spica
+  (test 19.1825 213.9154		;Arcturus
+	-11.1614 201.2983		;Spica
+	32.7930 0.00005)
+
+  ;; Front to nearly the back.
+  (test 0. 0.
+	0. 135.
+	135. 0.0000000000005)
+
+  ;; Across the N pole.
+  (test 80. 20.
+	80. 200.
+	20. 0.0000000000005)
+
+  ;; Top to nearly bottom on right.
+  (test 90. 90.
+	-85. 90.
+	175. 0.0000000000005))
+
 (define (2d-translate! p dx dy q)
   (let ((x (x p))
 	(y (y p)))
@@ -119,4 +191,9 @@ USA.
     (let ((x (x p))
 	  (y (y p)))
       (flo:vector-set! q 0 (flo:+ (flo:* x cos) (flo:* y (flo:negate sin))))
-      (flo:vector-set! q 1 (flo:+ (flo:* x sin) (flo:* y cos))))))
\ No newline at end of file
+      (flo:vector-set! q 1 (flo:+ (flo:* x sin) (flo:* y cos))))))
+
+(define (flo:~= n1 n2 epsilon)
+  ;; aka approximately-=
+  (and (flo:flonum? n1) (flo:flonum? n2) (flo:flonum? epsilon)
+       (flo:< (flo:abs (flo:- n2 n1)) epsilon)))
\ No newline at end of file
diff --git a/src/planetarium/mit.pkg b/src/planetarium/mit.pkg
index 56dcb8918..ad91156b5 100644
--- a/src/planetarium/mit.pkg
+++ b/src/planetarium/mit.pkg
@@ -98,9 +98,9 @@ USA.
 	  ;; and thus do not have to be imported, but are listed here
 	  ;; anyway for completeness (and analysis someday).
 	  guarantee-integer integer-divide round->exact truncate->exact
-	  ->flonum flo:< flo:<= flo:+ flo:- flo:* flo:/
+	  ->flonum flo:= flo:< flo:<= flo:+ flo:- flo:* flo:/
 	  flo:negative? flo:negate flo:truncate
-	  flo:sin flo:cos flo:atan2
+	  flo:sin flo:cos flo:atan2 flo:sqrt flo:atan
 	  flo:vector-cons flo:vector-ref flo:vector-set!
 
 	  create-thread detach-thread
diff --git a/src/planetarium/time.scm b/src/planetarium/time.scm
index 6cf9129ec..cfeb3c5dc 100644
--- a/src/planetarium/time.scm
+++ b/src/planetarium/time.scm
@@ -56,11 +56,6 @@ USA.
 	  B
 	  -1524.5)))))
 
-(define (flo:~= n1 n2 epsilon)
-  ;; aka approximately-=
-  (and (flo:flonum? n1) (flo:flonum? n2) (flo:flonum? epsilon)
-       (flo:< (flo:abs (flo:- n2 n1)) epsilon)))
-
 (define (test-julian-day)
   (for-each
     (lambda (item)
-- 
2.25.1