Merge pull request #4 from owlpinetech/feature/projection-update

Use flatsphere for more solid projection math

go.mod

@@ -3,3 +3,5 @@ module github.com/owlpinetech/healpix
 go 1.20
 
 require golang.org/x/exp v0.0.0-20230905200255-921286631fa9
+
+require github.com/owlpinetech/flatsphere v0.0.1 // indirect

go.sum

@@ -1,2 +1,4 @@
+github.com/owlpinetech/flatsphere v0.0.1 h1:rSkNlKx/YNzw6QyZcQhDGuLoN2RJugxDoyeTRThHD7E=
+github.com/owlpinetech/flatsphere v0.0.1/go.mod h1:Wc+2RhBlT09wCtZ1+Jc4CBYfUBnxs2vwSMaQRDFFWF8=
 golang.org/x/exp v0.0.0-20230905200255-921286631fa9 h1:GoHiUyI/Tp2nVkLI2mCxVkOjsbSXD66ic0XW0js0R9g=
 golang.org/x/exp v0.0.0-20230905200255-921286631fa9/go.mod h1:S2oDrQGGwySpoQPVqRShND87VCbxmc6bL1Yd2oYrm6k=

where.go

@@ -1,9 +1,9 @@
 package healpix
 
 import (
-	"fmt"
 	"math"
 
+	"github.com/owlpinetech/flatsphere"
 	"github.com/owlpinetech/healpix/internal/intmath"
 )
 
@@ -477,25 +477,8 @@ func (p ProjectionCoordinate) ToProjectionCoordinate(hp Healpix) ProjectionCoord
 }
 
 func (p ProjectionCoordinate) ToSphereCoordinate(hp Healpix) SphereCoordinate {
-	absY := math.Abs(p.y)
-	if absY >= math.Pi/2 {
-		panic(fmt.Sprintf("healpix: domain error in projection coordinate y dimension, %v too big", p.y))
-	}
-
-	if absY <= math.Pi/4 {
-		// equatorial region
-		z := (8 / (3 * math.Pi)) * p.y
-		colat := math.Acos(z)
-		return SphereCoordinate{math.Pi/2 - colat, colat, p.x}
-	} else {
-		// polar region
-		tt := math.Mod(p.x, math.Pi/2)
-		lng := p.x - ((absY-math.Pi/4)/(absY-math.Pi/2))*(tt-math.Pi/4)
-		zz := 2 - 4*absY/math.Pi
-		z := (1 - 1.0/3.0*(zz*zz)) * (p.y / absY)
-		colat := math.Acos(z)
-		return SphereCoordinate{math.Pi/2 - colat, colat, lng}
-	}
+	lat, lon := flatsphere.NewHEALPixStandard().Inverse(p.x, p.y)
+	return NewLatLonCoordinate(lat, lon)
 }
 
 func (p ProjectionCoordinate) PixelId(hp Healpix, scheme HealpixScheme) int {
@@ -560,21 +543,8 @@ func (p SphereCoordinate) ToRingCoordinate(hp Healpix) RingCoordinate {
 }
 
 func (p SphereCoordinate) ToProjectionCoordinate(hp Healpix) ProjectionCoordinate {
-	z := math.Cos(p.colatitude)
-	if math.Abs(z) <= 2.0/3.0 {
-		// equatiorial region
-		return ProjectionCoordinate{p.longitude, 3 * (math.Pi / 8) * z}
-	} else {
-		// polar region
-		facetX := math.Mod(p.longitude, math.Pi/2)
-		sigma := 2 - math.Sqrt(3*(1-math.Abs(z)))
-		if z < 0 {
-			sigma = -sigma
-		}
-		y := (math.Pi / 4) * sigma
-		x := p.longitude - (math.Abs(sigma)-1)*(facetX-math.Pi/4)
-		return ProjectionCoordinate{x, y}
-	}
+	x, y := flatsphere.NewHEALPixStandard().Project(p.Latitude(), p.Longitude())
+	return NewProjectionCoordinate(x, y)
 }
 
 func (p SphereCoordinate) ToSphereCoordinate(hp Healpix) SphereCoordinate {

where_test.go

@@ -228,37 +228,6 @@ func TestPositionToNestPixel(t *testing.T) {
 	}
 }
 
-func TestPositionToProjection(t *testing.T) {
-	testCases := []struct {
-		name       string
-		colatitude float64
-		longitude  float64
-		xproj      float64
-		yproj      float64
-	}{
-		{"lat/lng Pi/2,0 = x/y 0,0", math.Pi / 2, 0, 0, 0},
-		{"lat/lng Pi/2,2Pi = x/y 2Pi,0", math.Pi / 2, 2 * math.Pi, 2 * math.Pi, 0},
-		{"lat/lng arccos(1/3),Pi = x/y Pi,Pi/8", math.Acos(1.0 / 3.0), math.Pi, math.Pi, math.Pi / 8},
-		{"lat/lng Pi/4,Pi = x/y Pi,Pi/4", math.Acos(2.0 / 3.0), math.Pi, math.Pi, math.Pi / 4},
-	}
-
-	for _, tc := range testCases {
-		hp := HealpixOrder(0)
-		t.Run(tc.name, func(t *testing.T) {
-			pos := SphereCoordinate{math.Pi/2 - tc.colatitude, tc.colatitude, tc.longitude}
-			proj := ProjectionCoordinate{tc.xproj, tc.yproj}
-			rProj := pos.ToProjectionCoordinate(hp)
-			rPos := proj.ToSphereCoordinate(hp)
-			if !withinTolerance(rPos.colatitude, tc.colatitude, 0.000000001) || !withinTolerance(rPos.longitude, tc.longitude, 0.000000001) {
-				t.Errorf("Projection to position expected %v,%v but got %v,%v instead", tc.colatitude, tc.longitude, rPos.colatitude, rPos.longitude)
-			}
-			if !withinTolerance(rProj.x, tc.xproj, 0.000000001) || !withinTolerance(rProj.y, tc.yproj, 0.000000001) {
-				t.Errorf("Position to projection expected %v,%v but got %v,%v instead", tc.xproj, tc.yproj, rProj.x, rProj.y)
-			}
-		})
-	}
-}
-
 func TestProjectionToFacePixel(t *testing.T) {
 	testCases := []struct {
 		name  string
@@ -305,33 +274,6 @@ func TestProjectionToFacePixel(t *testing.T) {
 	}
 }
 
-func TestSphereToProjectionInvertible(t *testing.T) {
-	testCases := []struct {
-		name       string
-		colatitude float64
-		longitude  float64
-	}{
-		{"colat/lng Pi/2,0", math.Pi / 2, 0},
-		{"colat/lng Pi/4,Pi", math.Pi / 4, math.Pi},
-		{"colat/lng Pi/3,3Pi/8", math.Pi / 3, 3 * math.Pi / 8},
-		{"colat/lng Pi/6,7Pi/4", math.Pi / 6, 7 * math.Pi / 4},
-		{"colat/lng 1,1", 1, 1},
-		{"colat/lng 3.14,6.28", 3.14, 6.28},
-	}
-
-	for _, tc := range testCases {
-		// the healpix argument is not actually used in sphere-projection conversions
-		hp := HealpixOrder(0)
-		t.Run(tc.name, func(t *testing.T) {
-			pos := SphereCoordinate{math.Pi/2 - tc.colatitude, tc.colatitude, tc.longitude}
-			rPos := pos.ToProjectionCoordinate(hp).ToSphereCoordinate(hp)
-			if !withinTolerance(rPos.colatitude, tc.colatitude, 0.000000001) || !withinTolerance(rPos.longitude, tc.longitude, 0.000000001) {
-				t.Errorf("Position inverted expected %v,%v but got %v,%v instead", tc.colatitude, tc.longitude, rPos.colatitude, rPos.longitude)
-			}
-		})
-	}
-}
-
 func TestRingCoordToFacePixelInvertible(t *testing.T) {
 	testCases := []struct {
 		name        string