Tactical commit before a change I'm uncertain of

src/yyy_data/core.clj

@@ -118,6 +118,7 @@
   (geopoint [this datum] (vector3d-to-geopoint this datum))
   Transformable
   (apply-transform [this transform]
+                   ;; tests say this works
     (println (str "(apply-transform " this " " transform ")"))
                    (let
                      [s1 (+' (/ (:s transform) 1e6) 1)           ;; scale
@@ -151,6 +152,7 @@
   expects a map `t` whose values are numbers, and returns a map which has
   for each key in `t` a number which is the inverse of the value of the
   same key in `t`."
+  ;; tests say this works
   [t]
   (reduce
     merge
@@ -237,7 +239,7 @@
   (grid-y [location] n)
   (latitude [location] (latitude (geopoint location :WGS84)))
   (longitude [location] (longitude (geopoint location :WGS84)))
-  (vector3d [location])
+  (vector3d [location] (geopoint-to-vector3d (osgrid-to-geopoint location :WGS84)))
   (geopoint [location datum]
             (osgrid-to-geopoint location datum))
   (osgrid [location datum] location))
@@ -246,6 +248,7 @@
 (defn raw-vector3d-to-geopoint
   ;; My initial design decision to allow datums to be passed as keywords or as maps
   ;; was WRONG: a datum SHALL BE passed as a key
+  ;; Tests say this works
   ([v datum]
    (vector3d-to-geopoint (:x v) (:y v) (:z v) datum))
   ([x y z d]
@@ -281,21 +284,22 @@
    (geopoint-to-osgrid (:latitude gp) (:longitude gp) (:datum gp) datum))
   ([latitude longitude from-datum to-datum]
    ;; for bizarrely over-complicated trigonometry, look no further.
+   ;; This doesn't work. But, to be brutally honest, I don't need it to.
    (let [point (GeoPoint. latitude longitude to-datum)
          φ     (radians latitude)
          λ     (radians longitude)
-         a     6377563.396       ;; Airy 1830 major & minor semi-axes
+         a     6377563.396       ;; Airy 1830 major semi-axis
-         b     6356256.909
+         b     6356256.909       ;; Airy 1830 minor semi-axis
          F0    0.9996012717      ;; OS Grid scale factor on central meridian
          φ0    (radians 49)      ;; OS Grid true origin latitude
          λ0    (radians -2)      ;; OS Grid true origin longitude
          Δφ    (-' φ φ0)         ;; latitude offset from origin
          Δλ    (-' λ λ0)         ;; longitude offset from origin
-         Δλ²   (reduce * 1 (repeat 2 Δλ))
+         Δλ²   (expt Δλ 2)
-         Δλ³   (reduce * 1 (repeat 3 Δλ))
+         Δλ³   (expt Δλ 3)
-         Δλ⁴   (reduce * 1 (repeat 4 Δλ))
+         Δλ⁴   (expt Δλ 4)
-         Δλ⁵   (reduce * 1 (repeat 5 Δλ))
+         Δλ⁵   (expt Δλ 5)
-         Δλ⁶   (reduce * 1 (repeat 6 Δλ))
+         Δλ⁶   (expt Δλ 6)
          N0    -100000           ;; northing of true origin, metres
          E0    400000            ;; easting of true origin, metres
          e²    (-' 1 (/ (*' b b) (*' a a)))  ;; eccentricity squared
@@ -349,9 +353,10 @@
 
 
 (defn raw-geopoint-to-vector3d
+  ;; Tests say this no longer works
   ([gp]
    (println (str "(geopoint-to-vector3d " geopoint ")"))
-   (geopoint-to-vector3d (:latitude gp) (:longitude gp) (datums (:datum gp))))
+   (geopoint-to-vector3d (:latitude gp) (:longitude gp) (datum gp)))
   ([latitude longitude datum]
    (let [φ (radians latitude)
          λ (radians longitude)
@@ -374,9 +379,12 @@
 
 
 (defn raw-osgrid-to-geopoint
+  ([osgrid]
+   (osgrid-to-geopoint osgrid :WGS84))
   ([osgrid datum]
    (osgrid-to-geopoint (:e osgrid) (:n osgrid) datum))
   ([E N datum]
+   ;; this currently doesn't work
    (let
      [a     6377563.396         ;; Airy 1830 major semi-axis
       b     6356256.909         ;; Airy 1830 minor semi-axis
@@ -385,19 +393,25 @@
       λ0    (radians -2)        ;; national grid true origin longitude
       N0    -100000             ;; northing of true origin, metres
       E0    400000              ;; easting of true origin, metres
+      ;; var e2 = 1 - (b*b)/(a*a);
       e²    (-' 1
                (/
                  (*' b b)
                  (*' a a)))      ;; eccentricity squared
+      ;; var n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n;
       n     (/ (-' a b) (+' a b))
-      n²    (*' n n)
+      n²    (expt n 2)
-      n³    (*' n n n)
+      n³    (expt n 3)
       [M φ] (loop  [φ φ0 M 0]
               (let
+                ;; φ = (N-N0-M)/(a*F0) + φ;
                 [φ₁ (+' φ (/ (-' N N0 M) (*' a F0)))
+                 ;; (φ-φ0)
                  Δφ (-' φ φ0)
+                 ;; var Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (φ-φ0);
                  Ma (*' (+' 1 n (*' 5/4 n²) (*' 5/4 n³)) Δφ)
-                 Mb (*' (+' (*' n 3) (*' n² 3) (*' n³ 21/8)) (sin Δφ) (cos Δφ))
+                 ;; var Mb = (3*n + 3*n*n + (21/8)*n3) * Math.sin(φ-φ0) * Math.cos(φ+φ0);
+                 Mb (*' (+' (*' n 3) (*' n² 3) (*' 21/8 n³)) (sin Δφ) (cos (+' φ φ0)))
                  Mc (*' (+' (*' n² 15/8) (*' n³ 15/8)) (sin (*' 2 Δφ)) (cos (*' 2 Δφ)))
                  Md (*' 35/24 n³ (sin (*' 3 Δφ)) (cos (*' 3 Δφ)))
                  M₁ (*' b F0 (+' (-' 0 Ma) Mb Mc (-' 0 Md)))]

test/yyy_data/core_test.clj

@@ -53,13 +53,58 @@
     (let
       [hazelfield (yyy_data.core.OsGrid. 277656 549165)
        latlon (geopoint hazelfield :WGS84)
-       hazelfield' (osgrid latlon :WGS84)]
+       actual (osgrid latlon :WGS84)]
       (is
-        (< (abs (- (:e hazelfield) (:e hazelfield'))) 0.001)
+        (< (abs (- (:e hazelfield) (:e actual))) 0.001)
           (str "Components should be close to one another: easting: "
-               (:e hazelfield) ", " (:e hazelfield')))
+               (:e hazelfield) ", " (:e actual)))
       (is
-        (< (abs (- (:n hazelfield) (:n hazelfield'))) 0.001)
+        (< (abs (- (:n hazelfield) (:n actual))) 0.001)
           (str "Components should be close to one another: northing: "
-               (:n hazelfield) ", " (:n hazelfield'))))))
+               (:n hazelfield) ", " (:n actual))))))
 
+
+(deftest geopoint-to-vector-round-trip-tests
+  (testing "round trip geopoint->vector->geopoint"
+    (let
+      [hazelfield (yyy_data.core.GeoPoint. 54.822218 -3.906009 :WGS84)
+       v (vector3d hazelfield)
+       actual (geopoint v :WGS84)]
+      (is
+        (< (abs (- (:latitude actual) 54.822218)) 0.001)
+        (str "Latitude should be 54.822219: "  (:latitude actual)))
+      (is
+        (< (abs (- (:longitude actual) -3.906009)) 0.001)
+        (str "Longitude should be -3.906009: "  (:longitude actual))))))
+
+
+(deftest osgrid-to-vector-round-trip-tests
+  (testing "round trip osgrid->vector->osgrid"
+    (let
+      [hazelfield (yyy_data.core.OsGrid. 277656 549165)
+       v (vector3d hazelfield)
+       actual (osgrid v :WGS84)]
+      (is
+        (< (abs (- (:e hazelfield) (:e actual))) 0.001)
+        (str "Components should be close to one another: easting: "
+             (:e hazelfield) ", " (:e actual)))
+      (is
+        (< (abs (- (:n hazelfield) (:n actual))) 0.001)
+        (str "Components should be close to one another: northing: "
+             (:n hazelfield) ", " (:n actual))))))
+
+
+;; Currently blows up horribly (mutual recursion/stack). But I currently
+;; don't need it to work.
+;; (deftest geopoint-change-datum-round-trip-wgs84-tests
+;;   (testing "round trip geopoint/wgs84->geopoint/osgb36->geopoint/wgs84"
+;;     (let
+;;       [hazelfield (yyy_data.core.GeoPoint. 54.822218 -3.906009 :WGS84)
+;;        gp (geopoint hazelfield :OSGB36)
+;;        actual (geopoint gp :WGS84)]
+;;       (is
+;;         (< (abs (- (:latitude actual) 54.822218)) 0.001)
+;;         (str "Latitude should be 54.822219: "  (:latitude actual)))
+;;       (is
+;;         (< (abs (- (:longitude actual) -3.906009)) 0.001)
+;;         (str "Longitude should be -3.906009: "  (:longitude actual))))))