1 // Copyright 2015 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s1 16 17 import ( 18 "math" 19 ) 20 21 // ChordAngle represents the angle subtended by a chord (i.e., the straight 22 // line segment connecting two points on the sphere). Its representation 23 // makes it very efficient for computing and comparing distances, but unlike 24 // Angle it is only capable of representing angles between 0 and π radians. 25 // Generally, ChordAngle should only be used in loops where many angles need 26 // to be calculated and compared. Otherwise it is simpler to use Angle. 27 // 28 // ChordAngle loses some accuracy as the angle approaches π radians. 29 // There are several different ways to measure this error, including the 30 // representational error (i.e., how accurately ChordAngle can represent 31 // angles near π radians), the conversion error (i.e., how much precision is 32 // lost when an Angle is converted to an ChordAngle), and the measurement 33 // error (i.e., how accurate the ChordAngle(a, b) constructor is when the 34 // points A and B are separated by angles close to π radians). All of these 35 // errors differ by a small constant factor. 36 // 37 // For the measurement error (which is the largest of these errors and also 38 // the most important in practice), let the angle between A and B be (π - x) 39 // radians, i.e. A and B are within "x" radians of being antipodal. The 40 // corresponding chord length is 41 // 42 // r = 2 * sin((π - x) / 2) = 2 * cos(x / 2) 43 // 44 // For values of x not close to π the relative error in the squared chord 45 // length is at most 4.5 * dblEpsilon (see MaxPointError below). 46 // The relative error in "r" is thus at most 2.25 * dblEpsilon ~= 5e-16. To 47 // convert this error into an equivalent angle, we have 48 // 49 // |dr / dx| = sin(x / 2) 50 // 51 // and therefore 52 // 53 // |dx| = dr / sin(x / 2) 54 // = 5e-16 * (2 * cos(x / 2)) / sin(x / 2) 55 // = 1e-15 / tan(x / 2) 56 // 57 // The maximum error is attained when 58 // 59 // x = |dx| 60 // = 1e-15 / tan(x / 2) 61 // ~= 1e-15 / (x / 2) 62 // ~= sqrt(2e-15) 63 // 64 // In summary, the measurement error for an angle (π - x) is at most 65 // 66 // dx = min(1e-15 / tan(x / 2), sqrt(2e-15)) 67 // (~= min(2e-15 / x, sqrt(2e-15)) when x is small) 68 // 69 // On the Earth's surface (assuming a radius of 6371km), this corresponds to 70 // the following worst-case measurement errors: 71 // 72 // Accuracy: Unless antipodal to within: 73 // --------- --------------------------- 74 // 6.4 nanometers 10,000 km (90 degrees) 75 // 1 micrometer 81.2 kilometers 76 // 1 millimeter 81.2 meters 77 // 1 centimeter 8.12 meters 78 // 28.5 centimeters 28.5 centimeters 79 // 80 // The representational and conversion errors referred to earlier are somewhat 81 // smaller than this. For example, maximum distance between adjacent 82 // representable ChordAngle values is only 13.5 cm rather than 28.5 cm. To 83 // see this, observe that the closest representable value to r^2 = 4 is 84 // r^2 = 4 * (1 - dblEpsilon / 2). Thus r = 2 * (1 - dblEpsilon / 4) and 85 // the angle between these two representable values is 86 // 87 // x = 2 * acos(r / 2) 88 // = 2 * acos(1 - dblEpsilon / 4) 89 // ~= 2 * asin(sqrt(dblEpsilon / 2) 90 // ~= sqrt(2 * dblEpsilon) 91 // ~= 2.1e-8 92 // 93 // which is 13.5 cm on the Earth's surface. 94 // 95 // The worst case rounding error occurs when the value halfway between these 96 // two representable values is rounded up to 4. This halfway value is 97 // r^2 = (4 * (1 - dblEpsilon / 4)), thus r = 2 * (1 - dblEpsilon / 8) and 98 // the worst case rounding error is 99 // 100 // x = 2 * acos(r / 2) 101 // = 2 * acos(1 - dblEpsilon / 8) 102 // ~= 2 * asin(sqrt(dblEpsilon / 4) 103 // ~= sqrt(dblEpsilon) 104 // ~= 1.5e-8 105 // 106 // which is 9.5 cm on the Earth's surface. 107 type ChordAngle float64 108 109 const ( 110 // NegativeChordAngle represents a chord angle smaller than the zero angle. 111 // The only valid operations on a NegativeChordAngle are comparisons, 112 // Angle conversions, and Successor/Predecessor. 113 NegativeChordAngle = ChordAngle(-1) 114 115 // RightChordAngle represents a chord angle of 90 degrees (a "right angle"). 116 RightChordAngle = ChordAngle(2) 117 118 // StraightChordAngle represents a chord angle of 180 degrees (a "straight angle"). 119 // This is the maximum finite chord angle. 120 StraightChordAngle = ChordAngle(4) 121 122 // maxLength2 is the square of the maximum length allowed in a ChordAngle. 123 maxLength2 = 4.0 124 ) 125 126 // ChordAngleFromAngle returns a ChordAngle from the given Angle. 127 func ChordAngleFromAngle(a Angle) ChordAngle { 128 if a < 0 { 129 return NegativeChordAngle 130 } 131 if a.isInf() { 132 return InfChordAngle() 133 } 134 l := 2 * math.Sin(0.5*math.Min(math.Pi, a.Radians())) 135 return ChordAngle(l * l) 136 } 137 138 // ChordAngleFromSquaredLength returns a ChordAngle from the squared chord length. 139 // Note that the argument is automatically clamped to a maximum of 4 to 140 // handle possible roundoff errors. The argument must be non-negative. 141 func ChordAngleFromSquaredLength(length2 float64) ChordAngle { 142 if length2 > maxLength2 { 143 return StraightChordAngle 144 } 145 return ChordAngle(length2) 146 } 147 148 // Expanded returns a new ChordAngle that has been adjusted by the given error 149 // bound (which can be positive or negative). Error should be the value 150 // returned by either MaxPointError or MaxAngleError. For example: 151 // 152 // a := ChordAngleFromPoints(x, y) 153 // a1 := a.Expanded(a.MaxPointError()) 154 func (c ChordAngle) Expanded(e float64) ChordAngle { 155 // If the angle is special, don't change it. Otherwise clamp it to the valid range. 156 if c.isSpecial() { 157 return c 158 } 159 return ChordAngle(math.Max(0.0, math.Min(maxLength2, float64(c)+e))) 160 } 161 162 // Angle converts this ChordAngle to an Angle. 163 func (c ChordAngle) Angle() Angle { 164 if c < 0 { 165 return -1 * Radian 166 } 167 if c.IsInfinity() { 168 return InfAngle() 169 } 170 return Angle(2 * math.Asin(0.5*math.Sqrt(float64(c)))) 171 } 172 173 // InfChordAngle returns a chord angle larger than any finite chord angle. 174 // The only valid operations on an InfChordAngle are comparisons, Angle 175 // conversions, and Successor/Predecessor. 176 func InfChordAngle() ChordAngle { 177 return ChordAngle(math.Inf(1)) 178 } 179 180 // IsInfinity reports whether this ChordAngle is infinite. 181 func (c ChordAngle) IsInfinity() bool { 182 return math.IsInf(float64(c), 1) 183 } 184 185 // isSpecial reports whether this ChordAngle is one of the special cases. 186 func (c ChordAngle) isSpecial() bool { 187 return c < 0 || c.IsInfinity() 188 } 189 190 // isValid reports whether this ChordAngle is valid or not. 191 func (c ChordAngle) isValid() bool { 192 return (c >= 0 && c <= maxLength2) || c.isSpecial() 193 } 194 195 // Successor returns the smallest representable ChordAngle larger than this one. 196 // This can be used to convert a "<" comparison to a "<=" comparison. 197 // 198 // Note the following special cases: 199 // 200 // NegativeChordAngle.Successor == 0 201 // StraightChordAngle.Successor == InfChordAngle 202 // InfChordAngle.Successor == InfChordAngle 203 func (c ChordAngle) Successor() ChordAngle { 204 if c >= maxLength2 { 205 return InfChordAngle() 206 } 207 if c < 0 { 208 return 0 209 } 210 return ChordAngle(math.Nextafter(float64(c), 10.0)) 211 } 212 213 // Predecessor returns the largest representable ChordAngle less than this one. 214 // 215 // Note the following special cases: 216 // 217 // InfChordAngle.Predecessor == StraightChordAngle 218 // ChordAngle(0).Predecessor == NegativeChordAngle 219 // NegativeChordAngle.Predecessor == NegativeChordAngle 220 func (c ChordAngle) Predecessor() ChordAngle { 221 if c <= 0 { 222 return NegativeChordAngle 223 } 224 if c > maxLength2 { 225 return StraightChordAngle 226 } 227 228 return ChordAngle(math.Nextafter(float64(c), -10.0)) 229 } 230 231 // MaxPointError returns the maximum error size for a ChordAngle constructed 232 // from 2 Points x and y, assuming that x and y are normalized to within the 233 // bounds guaranteed by s2.Point.Normalize. The error is defined with respect to 234 // the true distance after the points are projected to lie exactly on the sphere. 235 func (c ChordAngle) MaxPointError() float64 { 236 // There is a relative error of (2.5*dblEpsilon) when computing the squared 237 // distance, plus a relative error of 2 * dblEpsilon, plus an absolute error 238 // of (16 * dblEpsilon**2) because the lengths of the input points may differ 239 // from 1 by up to (2*dblEpsilon) each. (This is the maximum error in Normalize). 240 return 4.5*dblEpsilon*float64(c) + 16*dblEpsilon*dblEpsilon 241 } 242 243 // MaxAngleError returns the maximum error for a ChordAngle constructed 244 // as an Angle distance. 245 func (c ChordAngle) MaxAngleError() float64 { 246 return dblEpsilon * float64(c) 247 } 248 249 // Add adds the other ChordAngle to this one and returns the resulting value. 250 // This method assumes the ChordAngles are not special. 251 func (c ChordAngle) Add(other ChordAngle) ChordAngle { 252 // Note that this method (and Sub) is much more efficient than converting 253 // the ChordAngle to an Angle and adding those and converting back. It 254 // requires only one square root plus a few additions and multiplications. 255 256 // Optimization for the common case where b is an error tolerance 257 // parameter that happens to be set to zero. 258 if other == 0 { 259 return c 260 } 261 262 // Clamp the angle sum to at most 180 degrees. 263 if c+other >= maxLength2 { 264 return StraightChordAngle 265 } 266 267 // Let a and b be the (non-squared) chord lengths, and let c = a+b. 268 // Let A, B, and C be the corresponding half-angles (a = 2*sin(A), etc). 269 // Then the formula below can be derived from c = 2 * sin(A+B) and the 270 // relationships sin(A+B) = sin(A)*cos(B) + sin(B)*cos(A) 271 // cos(X) = sqrt(1 - sin^2(X)) 272 x := float64(c * (1 - 0.25*other)) 273 y := float64(other * (1 - 0.25*c)) 274 return ChordAngle(math.Min(maxLength2, x+y+2*math.Sqrt(x*y))) 275 } 276 277 // Sub subtracts the other ChordAngle from this one and returns the resulting 278 // value. This method assumes the ChordAngles are not special. 279 func (c ChordAngle) Sub(other ChordAngle) ChordAngle { 280 if other == 0 { 281 return c 282 } 283 if c <= other { 284 return 0 285 } 286 x := float64(c * (1 - 0.25*other)) 287 y := float64(other * (1 - 0.25*c)) 288 return ChordAngle(math.Max(0.0, x+y-2*math.Sqrt(x*y))) 289 } 290 291 // Sin returns the sine of this chord angle. This method is more efficient 292 // than converting to Angle and performing the computation. 293 func (c ChordAngle) Sin() float64 { 294 return math.Sqrt(c.Sin2()) 295 } 296 297 // Sin2 returns the square of the sine of this chord angle. 298 // It is more efficient than Sin. 299 func (c ChordAngle) Sin2() float64 { 300 // Let a be the (non-squared) chord length, and let A be the corresponding 301 // half-angle (a = 2*sin(A)). The formula below can be derived from: 302 // sin(2*A) = 2 * sin(A) * cos(A) 303 // cos^2(A) = 1 - sin^2(A) 304 // This is much faster than converting to an angle and computing its sine. 305 return float64(c * (1 - 0.25*c)) 306 } 307 308 // Cos returns the cosine of this chord angle. This method is more efficient 309 // than converting to Angle and performing the computation. 310 func (c ChordAngle) Cos() float64 { 311 // cos(2*A) = cos^2(A) - sin^2(A) = 1 - 2*sin^2(A) 312 return float64(1 - 0.5*c) 313 } 314 315 // Tan returns the tangent of this chord angle. 316 func (c ChordAngle) Tan() float64 { 317 return c.Sin() / c.Cos() 318 } 319 320 // TODO(roberts): Differences from C++: 321 // Helpers to/from E5/E6/E7 322 // Helpers to/from degrees and radians directly. 323 // FastUpperBoundFrom(angle Angle) 324