@@ -77,4 +77,121 @@ import GeoInterface as GI
77
77
@test_throws AssertionError GeographicFromUnitSphere ()([1.0 , 0.0 ])
78
78
@test_throws AssertionError GeographicFromUnitSphere ()([1.0 , 0.0 , 0.0 , 0.0 ])
79
79
end
80
+ end
81
+
82
+ @testset " Spherical caps" begin
83
+ @testset " Construction" begin
84
+ # Test construction from UnitSphericalPoint
85
+ point = UnitSphericalPoint (1.0 , 0.0 , 0.0 )
86
+ cap = SphericalCap (point, π/ 4 )
87
+ @test cap. point == point
88
+ @test cap. radius == π/ 4
89
+
90
+ # Test construction from geographic point
91
+ geo_point = GI. Point (45 , 45 )
92
+ cap = SphericalCap (geo_point, π/ 4 )
93
+ @test cap. point isa UnitSphericalPoint{Float64}
94
+ @test cap. radius == π/ 4
95
+
96
+ # Test construction from tuple
97
+ cap = SphericalCap ((45 , 45 ), π/ 4 )
98
+ @test cap. point isa UnitSphericalPoint{Float64}
99
+ @test cap. radius == π/ 4
100
+ end
101
+
102
+ @testset " Intersection and containment" begin
103
+ # Create two caps that intersect
104
+ cap1 = SphericalCap (UnitSphericalPoint (1.0 , 0.0 , 0.0 ), π/ 4 )
105
+ cap2 = SphericalCap (UnitSphericalPoint (1 /√ 2 , 1 /√ 2 , 0.0 ), π/ 4 )
106
+ @test UnitSpherical. _intersects (cap1, cap2)
107
+ @test UnitSpherical. _intersects (cap2, cap1)
108
+ @test ! UnitSpherical. _disjoint (cap1, cap2)
109
+
110
+ # Create two caps that don't intersect
111
+ cap3 = SphericalCap (UnitSphericalPoint (1.0 , 0.0 , 0.0 ), π/ 8 )
112
+ cap4 = SphericalCap (UnitSphericalPoint (0.0 , 0.0 , 1.0 ), π/ 8 )
113
+ @test ! UnitSpherical. _intersects (cap3, cap4)
114
+ @test UnitSpherical. _disjoint (cap3, cap4)
115
+
116
+ # Test containment
117
+ big_cap = SphericalCap (UnitSphericalPoint (1.0 , 0.0 , 0.0 ), π/ 2 )
118
+ small_cap = SphericalCap (UnitSphericalPoint (1 /√ 2 , 1 /√ 2 , 0.0 ), π/ 4 )
119
+ @test UnitSpherical. _contains (big_cap, small_cap)
120
+ @test ! UnitSpherical. _contains (small_cap, big_cap)
121
+ end
122
+
123
+ @testset " Circumcenter and circumradius" begin
124
+ # Test with an equilateral triangle on the equator
125
+ a = UnitSphericalPoint (1.0 , 0.0 , 0.0 )
126
+ b = UnitSphericalPoint (- 0.5 , √ 3 / 2 , 0.0 )
127
+ c = UnitSphericalPoint (- 0.5 , -√ 3 / 2 , 0.0 )
128
+ cap = SphericalCap (a, b, c)
129
+
130
+ # The circumcenter should be at the north pole
131
+ @test isapprox (cap. point[1 ], 0.0 , atol= 1e-10 )
132
+ @test isapprox (cap. point[2 ], 0.0 , atol= 1e-10 )
133
+ @test isapprox (cap. point[3 ], 1.0 , atol= 1e-10 )
134
+ # The radius should be π/2 (90 degrees)
135
+ @test isapprox (cap. radius, π/ 2 , atol= 1e-10 )
136
+
137
+ # Test with a triangle in the northern hemisphere
138
+ a = UnitSphericalPoint (1.0 , 0.0 , 0.0 )
139
+ b = UnitSphericalPoint (0.0 , 1.0 , 0.0 )
140
+ c = UnitSphericalPoint (0.0 , 0.0 , 1.0 )
141
+ cap = SphericalCap (a, b, c)
142
+
143
+ # The circumcenter should be at (1/√3, 1/√3, 1/√3)
144
+ @test isapprox (cap. point[1 ], 1 /√ 3 , atol= 1e-10 )
145
+ @test isapprox (cap. point[2 ], 1 /√ 3 , atol= 1e-10 )
146
+ @test isapprox (cap. point[3 ], 1 /√ 3 , atol= 1e-10 )
147
+ # The radius should be the angle between the center and any vertex
148
+ expected_radius = acos (1 /√ 3 )
149
+ @test isapprox (cap. radius, expected_radius, atol= 1e-10 )
150
+
151
+ # Test with nearly colinear points (small angle between them)
152
+ # Points near the equator with very small angular separation
153
+ ϵ = 1e-6 # Very small angle in radians
154
+ a = UnitSphericalPoint (1.0 , 0.0 , 0.0 )
155
+ b = UnitSphericalPoint (cos (ϵ), sin (ϵ), 0.0 )
156
+ c = UnitSphericalPoint (cos (2 ϵ), sin (2 ϵ), 0.0 )
157
+ cap = SphericalCap (a, b, c)
158
+
159
+ # The circumcenter should be near the north pole
160
+ @test isapprox (cap. point[1 ], 0.0 , atol= 1e-6 )
161
+ @test isapprox (cap. point[2 ], 0.0 , atol= 1e-6 )
162
+ @test isapprox (cap. point[3 ], 1.0 , atol= 1e-6 )
163
+ # The radius should be approximately π/2
164
+ @test isapprox (cap. radius, π/ 2 , atol= 1e-6 )
165
+
166
+ # # Test with nearly identical points
167
+ # # Points very close to each other in the northern hemisphere
168
+ # ϵ = 1e-8 # Extremely small angle
169
+ # base = UnitSphericalPoint(1/√3, 1/√3, 1/√3)
170
+ # a = base
171
+ # b = UnitSphericalPoint(cos(ϵ), sin(ϵ), 0.0) * (1/√3) # Small perturbation
172
+ # c = UnitSphericalPoint(cos(2ϵ), sin(2ϵ), 0.0) * (1/√3) # Another small perturbation
173
+ # cap = SphericalCap(a, b, c)
174
+
175
+ # # The circumcenter should be very close to the base point
176
+ # @test isapprox(cap.point[1], 1/√3, atol=1e-6)
177
+ # @test isapprox(cap.point[2], 1/√3, atol=1e-6)
178
+ # @test isapprox(cap.point[3], 1/√3, atol=1e-6)
179
+ # # The radius should be very small
180
+ # @test isapprox(cap.radius, ϵ, atol=1e-6)
181
+
182
+ # Test with points forming a very thin triangle
183
+ # Points near the north pole with small angular separation
184
+ ϵ = 1e-6
185
+ a = UnitSphericalPoint (0.0 , 0.0 , 1.0 ) # North pole
186
+ b = UnitSphericalPoint (sin (ϵ), 0.0 , cos (ϵ))
187
+ c = UnitSphericalPoint (0.0 , sin (ϵ), cos (ϵ))
188
+ cap = SphericalCap (a, b, c)
189
+
190
+ # The circumcenter should be very close to the north pole
191
+ @test isapprox (cap. point[1 ], 0.0 , atol= 1e-6 )
192
+ @test isapprox (cap. point[2 ], 0.0 , atol= 1e-6 )
193
+ @test isapprox (cap. point[3 ], 1.0 , atol= 1e-6 )
194
+ # The radius should be approximately ϵ
195
+ @test isapprox (cap. radius, ϵ, atol= 1e-6 )
196
+ end
80
197
end
0 commit comments