Fix bugs in projection functions (#99)

* Add test case to catch issue #97

* Fix bugs in mollweideproj and equiproj

* [skipci] Update CHANGELOG

CHANGELOG.md

@@ -4,6 +4,8 @@
 
 -   Add functions `readClFromFITS`, `writeClToFITS`, `cl2dl`, `dl2cl`, `almxfl`, `almxfl!`, `almExplicitIndex` [#91](https://github.com/ziotom78/Healpix.jl/pull/91), thanks to [LeeoBianchi](https://github.com/LeeoBianchi)
 
+-   Fix a bug in `mollweideproj` and `equiproj` [#97](https://github.com/ziotom78/Healpix.jl/issues/97)
+
 # Version 4.1.2
 
 -   Use double precision in `ang2pixRing` and `zphi2pixRing` [#94](https://github.com/ziotom78/Healpix.jl/pull/94)

src/projections.jl

@@ -110,8 +110,8 @@ on the projection plane (both are in the range [−1, 1]).
 function equiproj(lat, lon)
     # We use `rem2pi` because we need angles in the range [-π, +π]
     x, y = (
-        (rem2pi(lon + π, RoundNearest) - π) / π,
-        2 * (rem2pi(lat + π, RoundNearest) - π) / π,
+        rem2pi(lon, RoundNearest) / π,
+        2lat / π,
     )
     (true, x, y)
 end
@@ -131,6 +131,17 @@ function equiprojinv(x, y)
     (true, π / 2 * y, π * x)
 end
 
+function find_mollweide_theta(ϕ; threshold = 1e-7)
+    abs(abs(ϕ) - π/2) < threshold && return ϕ
+
+    θ = ϕ
+    while true
+        new_θ = θ - (2θ + sin(2θ) - π * sin(ϕ)) / (2 + 2cos(2θ))
+        (abs(new_θ - θ) < threshold) && return θ
+        θ = new_θ
+    end
+end
+
 """
     mollweideproj(lat, lon)
 
@@ -141,8 +152,9 @@ or not (false), and the two numbers are the x and y coordinates of the point
 on the projection plane (both are in the range [−1, 1]).
 """
 function mollweideproj(lat, lon)
-    sinlat, coslat = sincos(lat)
-    (true, sinlat, 2 / π * lon * coslat)
+    θ = find_mollweide_theta(lat)
+
+    (true, -1 / π * lon * cos(θ), sin(θ))
 end
 
 """

test/test_projections.jl

@@ -1,6 +1,38 @@
 m = Healpix.HealpixMap{Float64,Healpix.RingOrder}(1)
 m.pixels = 1.0:12.0
 
+let projections = [
+    (Healpix.equiproj, Healpix.equiprojinv, "Equirectangular"),
+    (Healpix.mollweideproj, Healpix.mollweideprojinv, "Mollweide"),
+    ]
+    @testset "$name projection" for (idx, name) in enumerate((x[3] for x in projections))
+        @testset "Point ($x_in, $y_in)" for (x_in, y_in) in [
+            (0.0, 0.0),
+            (0.1, 0.1),
+            (0.9, 0.9),
+            (0.1, 0.9),
+            (0.9, 0.1),
+            (0.5, 0.1),
+            (0.5, 0.9),
+            (0.1, 0.5),
+            (0.9, 0.5),
+            ]
+
+            projfn = projections[idx][1]
+            invprojfn = projections[idx][2]
+
+            visible_inv, lat, long = invprojfn(x_in, y_in)
+            (!visible_inv) && continue
+            visible, x, y = projfn(lat, long)
+
+            @test visible_inv == visible
+
+            @test x ≈ x_in
+            @test y ≈ y_in
+        end
+    end
+end
+
 # Do not run @test here, just check that the function can be ran
 bmp = Healpix.equirectangular(m, Dict(:width => 50))
 bmp = Healpix.mollweide(m, Dict(:width => 50))