Merge pull request #43 from mweastwood/devel

fix interpolation and implement UNSEEN()

NEWS.md

@@ -1,20 +1,29 @@
 # News
 
+## v0.2.4
+
+*2017-09-15*
+
+* Fix a bug in the Julia port of interpolation.
+* Add `UNSEEN()` as a sentinal value for blind or masked pixels in Healpix maps. This value is
+  consistent with `healpy` and `libhealpix_cxx`.
+
 ## v0.2.3
 
 *2017-08-26*
 
-* implement `interpolate` in Julia (fixes a bug for interpolation near phi=2pi and provides a 30% speed boost)
-* use `verify_angles` in `query_disc` and `interpolate`
-* update to Healpix 3.31
+* Implement `interpolate` in Julia (fixes a bug for interpolation near phi=2pi and provides a 30%
+  speed boost).
+* Use `verify_angles` in `query_disc` and `interpolate`.
+* Update to Healpix 3.31.
 
 ## v0.2.2
 
 *2017-08-11*
 
-* implement `query_disc` for getting a list of all pixels interior to a disc
-* widen the signature of most pixel functions to `Integer`
-* fix the Mac build (`cfitsio` was moved away from the `science` tap)
+* Implement `query_disc` for getting a list of all pixels interior to a disc.
+* Widen the signature of most pixel functions to `Integer`.
+* Fix the Mac build (`cfitsio` was moved away from the `science` tap).
 
 ## v0.2.1
 

docs/src/library.md

@@ -39,6 +39,7 @@ ang2pix
 pix2ang
 vec2pix
 pix2vec
+LibHealpix.UNSEEN
 LibHealpix.interpolate
 query_disc
 ```

src/map.jl

@@ -341,6 +341,13 @@ end
 @inline _pixels(map::HealpixMap) = map.pixels
 @inline _pixels(other) = other
 
+"""
+    LibHealpix.UNSEEN()
+
+Get the sentinal value indicating a blind or masked pixel.
+"""
+UNSEEN() = -1.6375e30
+
 function interpolate_cxx(map::HealpixMap{Float32}, θ, ϕ)
     # this function is used to test the Julia implementation
     θ′, ϕ′ = verify_angles(θ, ϕ)
@@ -464,9 +471,10 @@ function interpolate(map, θ, ϕ)
     for idx = 1:4
         val = map[pix[idx]]
         wval = wgt[idx]
-        # if !unseen
-        numerator   += wval * val
-        denominator += wval
+        if val != UNSEEN()
+            numerator   += wval * val
+            denominator += wval
+        end
     end
     numerator / denominator
 end

src/rings.jl

@@ -57,17 +57,17 @@ function ring_info2(nside, ring)
         θ = atan2(sinθ, cosθ)
         ringpix = 4northring
         shifted = true
-        startpix = 2northring*northring - 1
+        startpix = 2northring*(northring - 1)
     else
         θ = acos((2nside-northring)*fact1)
         ringpix = 4nside
         shifted = ((northring-nside) & 1) == 0
-        startpix = ncap + (northring-nside)*ringpix + 1
+        startpix = ncap + (northring-nside)*ringpix
     end
     if northring != ring # southern hemisphere
         θ = π - θ
-        startpix = npix - startpix - ringpix + 2
+        startpix = npix - startpix - ringpix
     end
-    startpix, ringpix, θ, shifted
+    startpix+1, ringpix, θ, shifted
 end
 

test/rings.jl

@@ -33,7 +33,7 @@
     end
 
     @testset "ring_info2" begin
-        nside = 2
+        nside = 4
         npix  = nside2npix( nside)
         nring = nside2nring(nside)
         for ring = 1:nring