func TestCapExpanded(t *testing.T) {
	cap50 := CapFromCenterAngle(xAxisPt, 50.0*s1.Degree)
	cap51 := CapFromCenterAngle(xAxisPt, 51.0*s1.Degree)

	if !empty.Expanded(s1.Angle(fullHeight)).IsEmpty() {
		t.Error("Expanding empty cap should return an empty cap")
	}
	if !full.Expanded(s1.Angle(fullHeight)).IsFull() {
		t.Error("Expanding a full cap should return an full cap")
	}

	if !cap50.Expanded(0).ApproxEqual(cap50) {
		t.Error("Expanding a cap by 0° should be equal to the original")
	}
	if !cap50.Expanded(1 * s1.Degree).ApproxEqual(cap51) {
		t.Error("Expanding 50° by 1° should equal the 51° cap")
	}

	if cap50.Expanded(129.99 * s1.Degree).IsFull() {
		t.Error("Expanding 50° by 129.99° should not give a full cap")
	}
	if !cap50.Expanded(130.01 * s1.Degree).IsFull() {
		t.Error("Expanding 50° by 130.01° should give a full cap")
	}
}

func TestCapContainsCell(t *testing.T) {
	faceRadius := math.Atan(math.Sqrt2)
	for face := 0; face < 6; face++ {
		// The cell consisting of the entire face.
		rootCell := CellFromCellID(CellIDFromFace(face))

		// A leaf cell at the midpoint of the v=1 edge.
		edgeCell := CellFromPoint(Point{faceUVToXYZ(face, 0, 1-eps)})

		// A leaf cell at the u=1, v=1 corner
		cornerCell := CellFromPoint(Point{faceUVToXYZ(face, 1-eps, 1-eps)})

		// Quick check for full and empty caps.
		if !full.ContainsCell(rootCell) {
			t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", full, rootCell, false, true)
		}

		// Check intersections with the bounding caps of the leaf cells that are adjacent to
		// cornerCell along the Hilbert curve.  Because this corner is at (u=1,v=1), the curve
		// stays locally within the same cube face.
		first := cornerCell.id.Advance(-3)
		last := cornerCell.id.Advance(4)
		for id := first; id < last; id = id.Next() {
			c := CellFromCellID(id).CapBound()
			if got, want := c.ContainsCell(cornerCell), id == cornerCell.id; got != want {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", c, cornerCell, got, want)
			}
		}

		for capFace := 0; capFace < 6; capFace++ {
			// A cap that barely contains all of capFace.
			center := unitNorm(capFace)
			covering := CapFromCenterAngle(center, s1.Angle(faceRadius+eps))
			if got, want := covering.ContainsCell(rootCell), capFace == face; got != want {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, rootCell, got, want)
			}
			if got, want := covering.ContainsCell(edgeCell), center.Vector.Dot(edgeCell.id.Point().Vector) > 0.1; got != want {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, edgeCell, got, want)
			}
			if got, want := covering.ContainsCell(edgeCell), covering.IntersectsCell(edgeCell); got != want {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, edgeCell, got, want)
			}
			if got, want := covering.ContainsCell(cornerCell), capFace == face; got != want {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", covering, cornerCell, got, want)
			}

			// A cap that barely intersects the edges of capFace.
			bulging := CapFromCenterAngle(center, s1.Angle(math.Pi/4+eps))
			if bulging.ContainsCell(rootCell) {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", bulging, rootCell, true, false)
			}
			if got, want := bulging.ContainsCell(edgeCell), capFace == face; got != want {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", bulging, edgeCell, got, want)
			}
			if bulging.ContainsCell(cornerCell) {
				t.Errorf("Cap(%v).ContainsCell(%v) = %t; want = %t", bulging, cornerCell, true, false)
			}
		}
	}
}

func TestCapContains(t *testing.T) {
	tests := []struct {
		c1, c2 Cap
		want   bool
	}{
		{empty, empty, true},
		{full, empty, true},
		{full, full, true},
		{empty, xAxis, false},
		{full, xAxis, true},
		{xAxis, full, false},
		{xAxis, xAxis, true},
		{xAxis, empty, true},
		{hemi, tiny, true},
		{hemi, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/4-epsilon)), true},
		{hemi, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/4+epsilon)), false},
		{concave, hemi, true},
		{concave, CapFromCenterHeight(Point{concave.center.Mul(-1.0)}, 0.1), false},
	}
	for _, test := range tests {
		if got := test.c1.Contains(test.c2); got != test.want {
			t.Errorf("%v.Contains(%v) = %t; want %t", test.c1, test.c2, got, test.want)
		}
	}
}

// Radius returns the cap's radius.
func (c Cap) Radius() s1.Angle {
	if c.IsEmpty() {
		return s1.Angle(emptyHeight)
	}

	// This could also be computed as acos(1 - height_), but the following
	// formula is much more accurate when the cap height is small. It
	// follows from the relationship h = 1 - cos(r) = 2 sin^2(r/2).
	return s1.Angle(2 * math.Asin(math.Sqrt(0.5*c.height)))
}

// samplePointFromRect returns a point chosen uniformly at random (with respect
// to area on the sphere) from the given rectangle.
func samplePointFromRect(rect Rect) Point {
	// First choose a latitude uniformly with respect to area on the sphere.
	sinLo := math.Sin(rect.Lat.Lo)
	sinHi := math.Sin(rect.Lat.Hi)
	lat := math.Asin(randomUniformFloat64(sinLo, sinHi))

	// Now choose longitude uniformly within the given range.
	lng := rect.Lng.Lo + randomFloat64()*rect.Lng.Length()

	return PointFromLatLng(LatLng{s1.Angle(lat), s1.Angle(lng)}.Normalized())
}

func TestCapAddPoint(t *testing.T) {
	tests := []struct {
		have Cap
		p    Point
		want Cap
	}{
		// Cap plus its center equals itself.
		{xAxis, xAxisPt, xAxis},
		{yAxis, yAxisPt, yAxis},

		// Cap plus opposite point equals full.
		{xAxis, PointFromCoords(-1, 0, 0), full},
		{yAxis, PointFromCoords(0, -1, 0), full},

		// Cap plus orthogonal axis equals half cap.
		{xAxis, PointFromCoords(0, 0, 1), CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/2.0))},
		{xAxis, PointFromCoords(0, 0, -1), CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/2.0))},

		// The 45 degree angled hemisphere plus some points.
		{
			hemi,
			PointFromCoords(0, 1, -1),
			CapFromCenterAngle(Point{PointFromCoords(1, 0, 1).Normalize()},
				s1.Angle(120.0)*s1.Degree),
		},
		{
			hemi,
			PointFromCoords(0, -1, -1),
			CapFromCenterAngle(Point{PointFromCoords(1, 0, 1).Normalize()},
				s1.Angle(120.0)*s1.Degree),
		},
		{
			// This angle between this point and the center is acos(-sqrt(2/3))
			hemi,
			PointFromCoords(-1, -1, -1),
			CapFromCenterAngle(Point{PointFromCoords(1, 0, 1).Normalize()},
				s1.Angle(2.5261129449194)),
		},
		{hemi, PointFromCoords(0, 1, 1), hemi},
		{hemi, PointFromCoords(1, 0, 0), hemi},
	}

	for _, test := range tests {
		got := test.have.AddPoint(test.p)
		if !got.ApproxEqual(test.want) {
			t.Errorf("%v.AddPoint(%v) = %v, want %v", test.have, test.p, got, test.want)
		}

		if !got.ContainsPoint(test.p) {
			t.Errorf("%v.AddPoint(%v) did not contain added point", test.have, test.p)
		}
	}
}

func rectFromDegrees(latLo, lngLo, latHi, lngHi float64) Rect {
	// Convenience method to construct a rectangle. This method is
	// intentionally *not* in the S2LatLngRect interface because the
	// argument order is ambiguous, but is fine for the test.
	return Rect{
		Lat: r1.Interval{
			Lo: (s1.Angle(latLo) * s1.Degree).Radians(),
			Hi: (s1.Angle(latHi) * s1.Degree).Radians(),
		},
		Lng: s1.IntervalFromEndpoints(
			(s1.Angle(lngLo) * s1.Degree).Radians(),
			(s1.Angle(lngHi) * s1.Degree).Radians(),
		),
	}
}

func TestCapAddCap(t *testing.T) {
	tests := []struct {
		have  Cap
		other Cap
		want  Cap
	}{
		// Identity cases.
		{empty, empty, empty},
		{full, full, full},

		// Anything plus empty equals itself.
		{full, empty, full},
		{empty, full, full},
		{xAxis, empty, xAxis},
		{empty, xAxis, xAxis},
		{yAxis, empty, yAxis},
		{empty, yAxis, yAxis},

		// Two halves make a whole.
		{xAxis, xComp, full},

		// Two zero-height orthogonal axis caps make a half-cap.
		{xAxis, yAxis, CapFromCenterAngle(xAxisPt, s1.Angle(math.Pi/2.0))},
	}

	for _, test := range tests {
		got := test.have.AddCap(test.other)
		if !got.ApproxEqual(test.want) {
			t.Errorf("%v.AddCap(%v) = %v, want %v", test.have, test.other, got, test.want)
		}
	}
}

// CapBound returns a cap that countains Rect.
func (r Rect) CapBound() Cap {
	// We consider two possible bounding caps, one whose axis passes
	// through the center of the lat-long rectangle and one whose axis
	// is the north or south pole.  We return the smaller of the two caps.

	if r.IsEmpty() {
		return EmptyCap()
	}

	var poleZ, poleAngle float64
	if r.Lat.Hi+r.Lat.Lo < 0 {
		// South pole axis yields smaller cap.
		poleZ = -1
		poleAngle = math.Pi/2 + r.Lat.Hi
	} else {
		poleZ = 1
		poleAngle = math.Pi/2 - r.Lat.Lo
	}
	poleCap := CapFromCenterAngle(PointFromCoords(0, 0, poleZ), s1.Angle(poleAngle)*s1.Radian)

	// For bounding rectangles that span 180 degrees or less in longitude, the
	// maximum cap size is achieved at one of the rectangle vertices.  For
	// rectangles that are larger than 180 degrees, we punt and always return a
	// bounding cap centered at one of the two poles.
	if math.Remainder(r.Lng.Hi-r.Lng.Lo, 2*math.Pi) >= 0 && r.Lng.Hi-r.Lng.Lo < 2*math.Pi {
		midCap := CapFromPoint(PointFromLatLng(r.Center())).AddPoint(PointFromLatLng(r.Lo())).AddPoint(PointFromLatLng(r.Hi()))
		if midCap.Height() < poleCap.Height() {
			return midCap
		}
	}
	return poleCap
}

// Distance returns the angle between two LatLngs.
func (ll LatLng) Distance(ll2 LatLng) s1.Angle {
	// Haversine formula, as used in C++ S2LatLng::GetDistance.
	lat1, lat2 := ll.Lat.Radians(), ll2.Lat.Radians()
	lng1, lng2 := ll.Lng.Radians(), ll2.Lng.Radians()
	dlat := math.Sin(0.5 * (lat2 - lat1))
	dlng := math.Sin(0.5 * (lng2 - lng1))
	x := dlat*dlat + dlng*dlng*math.Cos(lat1)*math.Cos(lat2)
	return s1.Angle(2*math.Atan2(math.Sqrt(x), math.Sqrt(math.Max(0, 1-x)))) * s1.Radian
}

// Normalized returns the normalized version of the LatLng,
// with Lat clamped to [-π/2,π/2] and Lng wrapped in [-π,π].
func (ll LatLng) Normalized() LatLng {
	lat := ll.Lat
	if lat > northPoleLat {
		lat = northPoleLat
	} else if lat < southPoleLat {
		lat = southPoleLat
	}
	lng := s1.Angle(math.Remainder(ll.Lng.Radians(), 2*math.Pi)) * s1.Radian
	return LatLng{lat, lng}
}

// Vertex returns the i-th vertex of the rectangle (i = 0,1,2,3) in CCW order
// (lower left, lower right, upper right, upper left).
func (r Rect) Vertex(i int) LatLng {
	var lat, lng float64

	switch i {
	case 0:
		lat = r.Lat.Lo
		lng = r.Lng.Lo
	case 1:
		lat = r.Lat.Lo
		lng = r.Lng.Hi
	case 2:
		lat = r.Lat.Hi
		lng = r.Lng.Hi
	case 3:
		lat = r.Lat.Hi
		lng = r.Lng.Lo
	}
	return LatLng{s1.Angle(lat) * s1.Radian, s1.Angle(lng) * s1.Radian}
}

// Interpolate returns the point X along the line segment AB whose distance from A
// is the given fraction "t" of the distance AB. Does NOT require that "t" be
// between 0 and 1. Note that all distances are measured on the surface of
// the sphere, so this is more complicated than just computing (1-t)*a + t*b
// and normalizing the result.
func Interpolate(t float64, a, b Point) Point {
	if t == 0 {
		return a
	}
	if t == 1 {
		return b
	}
	ab := a.Angle(b.Vector)
	return InterpolateAtDistance(s1.Angle(t)*ab, a, b)
}

func TestInterpolateOverLongEdge(t *testing.T) {
	lng := math.Pi - 1e-2
	a := Point{PointFromLatLng(LatLng{0, 0}).Normalize()}
	b := Point{PointFromLatLng(LatLng{0, s1.Angle(lng)}).Normalize()}

	for f := 0.4; f > 1e-15; f *= 0.1 {
		// Test that interpolation is accurate on a long edge (but not so long that
		// the definition of the edge itself becomes too unstable).
		want := Point{PointFromLatLng(LatLng{0, s1.Angle(f * lng)}).Normalize()}
		if got := Interpolate(f, a, b); !pointsApproxEquals(got, want, 3e-15) {
			t.Errorf("long edge Interpolate(%v, %v, %v) = %v, want %v", f, a, b, got, want)
		}

		// Test the remainder of the dist also matches.
		wantRem := Point{PointFromLatLng(LatLng{0, s1.Angle((1 - f) * lng)}).Normalize()}
		if got := Interpolate(1-f, a, b); !pointsApproxEquals(got, wantRem, 3e-15) {
			t.Errorf("long edge Interpolate(%v, %v, %v) = %v, want %v", 1-f, a, b, got, wantRem)
		}
	}
}

func TestRadiusToHeight(t *testing.T) {
	tests := []struct {
		got  s1.Angle
		want float64
	}{
		// Above/below boundary checks.
		{s1.Angle(-0.5), emptyHeight},
		{s1.Angle(0), 0},
		{s1.Angle(math.Pi), fullHeight},
		{s1.Angle(2 * math.Pi), fullHeight},
		// Degree tests.
		{-7.0 * s1.Degree, emptyHeight},
		{-0.0 * s1.Degree, 0},
		{0.0 * s1.Degree, 0},
		{12.0 * s1.Degree, 0.02185239926619},
		{30.0 * s1.Degree, 0.13397459621556},
		{45.0 * s1.Degree, 0.29289321881345},
		{90.0 * s1.Degree, 1.0},
		{179.99 * s1.Degree, 1.99999998476912},
		{180.0 * s1.Degree, fullHeight},
		{270.0 * s1.Degree, fullHeight},
		// Radians tests.
		{-1.0 * s1.Radian, emptyHeight},
		{-0.0 * s1.Radian, 0},
		{0.0 * s1.Radian, 0},
		{1.0 * s1.Radian, 0.45969769413186},
		{math.Pi / 2.0 * s1.Radian, 1.0},
		{2.0 * s1.Radian, 1.41614683654714},
		{3.0 * s1.Radian, 1.98999249660044},
		{math.Pi * s1.Radian, fullHeight},
		{4.0 * s1.Radian, fullHeight},
	}
	for _, test := range tests {
		// float64Eq comes from s2latlng_test.go
		if got := radiusToHeight(test.got); !float64Eq(got, test.want) {
			t.Errorf("radiusToHeight(%v) = %v; want %v", test.got, got, test.want)
		}
	}
}

func TestContainsPointConsistentWithS2CellIDFromPoint(t *testing.T) {
	// Construct many points that are nearly on a Cell edge, and verify that
	// CellFromCellID(cellIDFromPoint(p)).Contains(p) is always true.
	for iter := 0; iter < 1000; iter++ {
		cell := CellFromCellID(randomCellID())
		i1 := randomUniformInt(4)
		i2 := (i1 + 1) & 3
		v1 := cell.Vertex(i1)
		v2 := samplePointFromCap(CapFromCenterAngle(cell.Vertex(i2), s1.Angle(epsilon)))
		p := Interpolate(randomFloat64(), v1, v2)
		if !CellFromCellID(cellIDFromPoint(p)).ContainsPoint(p) {
			t.Errorf("For p=%v, CellFromCellID(cellIDFromPoint(p)).ContainsPoint(p) was false", p)
		}
	}
}

func TestRectCapBound(t *testing.T) {
	tests := []struct {
		r    Rect
		want Cap
	}{
		{ // Bounding cap at center is smaller.
			rectFromDegrees(-45, -45, 45, 45),
			CapFromCenterHeight(PointFromCoords(1, 0, 0), 0.5),
		},
		{ // Bounding cap at north pole is smaller.
			rectFromDegrees(88, -80, 89, 80),
			CapFromCenterAngle(PointFromCoords(0, 0, 1), s1.Angle(2)*s1.Degree),
		},
		{ // Longitude span > 180 degrees.
			rectFromDegrees(-30, -150, -10, 50),
			CapFromCenterAngle(PointFromCoords(0, 0, -1), s1.Angle(80)*s1.Degree),
		},
	}
	for _, test := range tests {
		if got := test.r.CapBound(); !test.want.ApproxEqual(got) {
			t.Errorf("%v.CapBound() was %v, want %v", test.r, got, test.want)
		}
	}
}

// perturbATowardsB returns a point that has been shifted some distance towards the
// second point based on a random number.
func perturbATowardsB(a, b Point) Point {
	choice := randomFloat64()
	if choice < 0.1 {
		return a
	}
	if choice < 0.3 {
		// Return a point that is exactly proportional to A and that still
		// satisfies IsUnitLength().
		for {
			b := Point{a.Mul(2 - a.Norm() + 5*(randomFloat64()-0.5)*dblEpsilon)}
			if !b.ApproxEqual(a) && b.IsUnit() {
				return b
			}
		}
	}
	if choice < 0.5 {
		// Return a point such that the distance squared to A will underflow.
		return InterpolateAtDistance(1e-300, a, b)
	}
	// Otherwise return a point whose distance from A is near dblEpsilon such
	// that the log of the pdf is uniformly distributed.
	distance := dblEpsilon * 1e-5 * math.Pow(1e6, randomFloat64())
	return InterpolateAtDistance(s1.Angle(distance), a, b)
}

func latitude(p Point) s1.Angle {
	return s1.Angle(math.Atan2(p.Z, math.Sqrt(p.X*p.X+p.Y*p.Y))) * s1.Radian
}

func TestExpandForSubregions(t *testing.T) {
	// Test the full and empty bounds.
	if !ExpandForSubregions(FullRect()).IsFull() {
		t.Errorf("Subregion Bound of full rect should be full")
	}
	if !ExpandForSubregions(EmptyRect()).IsEmpty() {
		t.Errorf("Subregion Bound of empty rect should be empty")
	}

	tests := []struct {
		xLat, xLng, yLat, yLng float64
		wantFull               bool
	}{
		// Cases where the bound does not straddle the equator (but almost does),
		// and spans nearly 180 degrees in longitude.
		{3e-16, 0, 1e-14, math.Pi, true},
		{9e-16, 0, 1e-14, math.Pi, false},
		{1e-16, 7e-16, 1e-14, math.Pi, true},
		{3e-16, 14e-16, 1e-14, math.Pi, false},
		{1e-100, 14e-16, 1e-14, math.Pi, true},
		{1e-100, 22e-16, 1e-14, math.Pi, false},
		// Cases where the bound spans at most 90 degrees in longitude, and almost
		// 180 degrees in latitude.  Note that DBL_EPSILON is about 2.22e-16, which
		// implies that the double-precision value just below Pi/2 can be written as
		// (math.Pi/2 - 2e-16).
		{-math.Pi / 2, -1e-15, math.Pi/2 - 7e-16, 0, true},
		{-math.Pi / 2, -1e-15, math.Pi/2 - 30e-16, 0, false},
		{-math.Pi/2 + 4e-16, 0, math.Pi/2 - 2e-16, 1e-7, true},
		{-math.Pi/2 + 30e-16, 0, math.Pi / 2, 1e-7, false},
		{-math.Pi/2 + 4e-16, 0, math.Pi/2 - 4e-16, math.Pi / 2, true},
		{-math.Pi / 2, 0, math.Pi/2 - 30e-16, math.Pi / 2, false},
		// Cases where the bound straddles the equator and spans more than 90
		// degrees in longitude.  These are the cases where the critical distance is
		// between a corner of the bound and the opposite longitudinal edge.  Unlike
		// the cases above, here the bound may contain nearly-antipodal points (to
		// within 3.055 * DBL_EPSILON) even though the latitude and longitude ranges
		// are both significantly less than (math.Pi - 3.055 * DBL_EPSILON).
		{-math.Pi / 2, 0, math.Pi/2 - 1e-8, math.Pi - 1e-7, true},
		{-math.Pi / 2, 0, math.Pi/2 - 1e-7, math.Pi - 1e-7, false},
		{-math.Pi/2 + 1e-12, -math.Pi + 1e-4, math.Pi / 2, 0, true},
		{-math.Pi/2 + 1e-11, -math.Pi + 1e-4, math.Pi / 2, 0, true},
	}

	for _, tc := range tests {
		in := RectFromLatLng(LatLng{s1.Angle(tc.xLat), s1.Angle(tc.xLng)})
		in = in.AddPoint(LatLng{s1.Angle(tc.yLat), s1.Angle(tc.yLng)})
		got := ExpandForSubregions(in)

		// Test that the bound is actually expanded.
		if !got.Contains(in) {
			t.Errorf("Subregion bound of (%f, %f, %f, %f) should contain original rect", tc.xLat, tc.xLng, tc.yLat, tc.yLng)
		}
		if in.Lat == validRectLatRange && in.Lat.ContainsInterval(got.Lat) {
			t.Errorf("Subregion bound of (%f, %f, %f, %f) shouldn't be contained by original rect", tc.xLat, tc.xLng, tc.yLat, tc.yLng)
		}

		// We check the various situations where the bound contains nearly-antipodal points. The tests are organized into pairs
		// where the two bounds are similar except that the first bound meets the nearly-antipodal criteria while the second does not.
		if got.IsFull() != tc.wantFull {
			t.Errorf("Subregion Bound of (%f, %f, %f, %f).IsFull should be %t", tc.xLat, tc.xLng, tc.yLat, tc.yLng, tc.wantFull)
		}
	}

	rectTests := []struct {
		xLat, xLng, yLat, yLng float64
		wantRect               Rect
	}{
		{1.5, -math.Pi / 2, 1.5, math.Pi/2 - 2e-16, Rect{r1.Interval{1.5, 1.5}, s1.FullInterval()}},
		{1.5, -math.Pi / 2, 1.5, math.Pi/2 - 7e-16, Rect{r1.Interval{1.5, 1.5}, s1.Interval{-math.Pi / 2, math.Pi/2 - 7e-16}}},
		// Check for cases where the bound is expanded to include one of the poles
		{-math.Pi/2 + 1e-15, 0, -math.Pi/2 + 1e-15, 0, Rect{r1.Interval{-math.Pi / 2, -math.Pi/2 + 1e-15}, s1.FullInterval()}},
		{math.Pi/2 - 1e-15, 0, math.Pi/2 - 1e-15, 0, Rect{r1.Interval{math.Pi/2 - 1e-15, math.Pi / 2}, s1.FullInterval()}},
	}

	for _, tc := range rectTests {
		// Now we test cases where the bound does not contain nearly-antipodal
		// points, but it does contain points that are approximately 180 degrees
		// apart in latitude.
		in := RectFromLatLng(LatLng{s1.Angle(tc.xLat), s1.Angle(tc.xLng)})
		in = in.AddPoint(LatLng{s1.Angle(tc.yLat), s1.Angle(tc.yLng)})
		got := ExpandForSubregions(in)
		if !rectsApproxEqual(got, tc.wantRect, rectErrorLat, rectErrorLng) {
			t.Errorf("Subregion Bound of (%f, %f, %f, %f) = (%v) should be %v", tc.xLat, tc.xLng, tc.yLat, tc.yLng, got, tc.wantRect)
		}
	}
}