planetarium: Add geodesic-distance, angular-separation,...

...and an angular-separation test procedure.

diff --git a/src/planetarium/geometry.scm b/src/planetarium/geometry.scm
index 45a9f90dcbbd6ad36d9b59c284bd971e5d04187e..9026b05c72354455fb374f61239c69185d65b412 100644 (file)
--- 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 56dcb8918560eb20dc37ca03f164040e976780b1..ad91156b547dd1fd0e9ca531f5d54f057e805893 100644 (file)
--- 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 6cf9129ecda68be613504f925cdf7aabd03c096d..cfeb3c5dc8a745d91414e034c2059807ff7f89f6 100644 (file)
--- 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)