1// Copyright (c) 2017 Couchbase, Inc. 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 15package geo 16 17import ( 18 "math" 19) 20 21var earthDiameterPerLatitude []float64 22var sinTab []float64 23var cosTab []float64 24var asinTab []float64 25var asinDer1DivF1Tab []float64 26var asinDer2DivF2Tab []float64 27var asinDer3DivF3Tab []float64 28var asinDer4DivF4Tab []float64 29 30const radiusTabsSize = (1 << 10) + 1 31const radiusDelta = (math.Pi / 2) / (radiusTabsSize - 1) 32const radiusIndexer = 1 / radiusDelta 33const sinCosTabsSize = (1 << 11) + 1 34const asinTabsSize = (1 << 13) + 1 35const oneDivF2 = 1 / 2.0 36const oneDivF3 = 1 / 6.0 37const oneDivF4 = 1 / 24.0 38 39// 1.57079632673412561417e+00 first 33 bits of pi/2 40var pio2Hi = math.Float64frombits(0x3FF921FB54400000) 41 42// 6.07710050650619224932e-11 pi/2 - PIO2_HI 43var pio2Lo = math.Float64frombits(0x3DD0B4611A626331) 44 45var asinPio2Hi = math.Float64frombits(0x3FF921FB54442D18) // 1.57079632679489655800e+00 46var asinPio2Lo = math.Float64frombits(0x3C91A62633145C07) // 6.12323399573676603587e-17 47var asinPs0 = math.Float64frombits(0x3fc5555555555555) // 1.66666666666666657415e-01 48var asinPs1 = math.Float64frombits(0xbfd4d61203eb6f7d) // -3.25565818622400915405e-01 49var asinPs2 = math.Float64frombits(0x3fc9c1550e884455) // 2.01212532134862925881e-01 50var asinPs3 = math.Float64frombits(0xbfa48228b5688f3b) // -4.00555345006794114027e-02 51var asinPs4 = math.Float64frombits(0x3f49efe07501b288) // 7.91534994289814532176e-04 52var asinPs5 = math.Float64frombits(0x3f023de10dfdf709) // 3.47933107596021167570e-05 53var asinQs1 = math.Float64frombits(0xc0033a271c8a2d4b) // -2.40339491173441421878e+00 54var asinQs2 = math.Float64frombits(0x40002ae59c598ac8) // 2.02094576023350569471e+00 55var asinQs3 = math.Float64frombits(0xbfe6066c1b8d0159) // -6.88283971605453293030e-01 56var asinQs4 = math.Float64frombits(0x3fb3b8c5b12e9282) // 7.70381505559019352791e-02 57 58var twoPiHi = 4 * pio2Hi 59var twoPiLo = 4 * pio2Lo 60var sinCosDeltaHi = twoPiHi/sinCosTabsSize - 1 61var sinCosDeltaLo = twoPiLo/sinCosTabsSize - 1 62var sinCosIndexer = 1 / (sinCosDeltaHi + sinCosDeltaLo) 63var sinCosMaxValueForIntModulo = ((math.MaxInt64 >> 9) / sinCosIndexer) * 0.99 64var asinMaxValueForTabs = math.Sin(73.0 * degreesToRadian) 65 66var asinDelta = asinMaxValueForTabs / (asinTabsSize - 1) 67var asinIndexer = 1 / asinDelta 68 69func init() { 70 // initializes the tables used for the sloppy math functions 71 72 // sin and cos 73 sinTab = make([]float64, sinCosTabsSize) 74 cosTab = make([]float64, sinCosTabsSize) 75 sinCosPiIndex := (sinCosTabsSize - 1) / 2 76 sinCosPiMul2Index := 2 * sinCosPiIndex 77 sinCosPiMul05Index := sinCosPiIndex / 2 78 sinCosPiMul15Index := 3 * sinCosPiIndex / 2 79 for i := 0; i < sinCosTabsSize; i++ { 80 // angle: in [0,2*PI]. 81 angle := float64(i)*sinCosDeltaHi + float64(i)*sinCosDeltaLo 82 sinAngle := math.Sin(angle) 83 cosAngle := math.Cos(angle) 84 // For indexes corresponding to null cosine or sine, we make sure the value is zero 85 // and not an epsilon. This allows for a much better accuracy for results close to zero. 86 if i == sinCosPiIndex { 87 sinAngle = 0.0 88 } else if i == sinCosPiMul2Index { 89 sinAngle = 0.0 90 } else if i == sinCosPiMul05Index { 91 sinAngle = 0.0 92 } else if i == sinCosPiMul15Index { 93 sinAngle = 0.0 94 } 95 sinTab[i] = sinAngle 96 cosTab[i] = cosAngle 97 } 98 99 // asin 100 asinTab = make([]float64, asinTabsSize) 101 asinDer1DivF1Tab = make([]float64, asinTabsSize) 102 asinDer2DivF2Tab = make([]float64, asinTabsSize) 103 asinDer3DivF3Tab = make([]float64, asinTabsSize) 104 asinDer4DivF4Tab = make([]float64, asinTabsSize) 105 for i := 0; i < asinTabsSize; i++ { 106 // x: in [0,ASIN_MAX_VALUE_FOR_TABS]. 107 x := float64(i) * asinDelta 108 asinTab[i] = math.Asin(x) 109 oneMinusXSqInv := 1.0 / (1 - x*x) 110 oneMinusXSqInv05 := math.Sqrt(oneMinusXSqInv) 111 oneMinusXSqInv15 := oneMinusXSqInv05 * oneMinusXSqInv 112 oneMinusXSqInv25 := oneMinusXSqInv15 * oneMinusXSqInv 113 oneMinusXSqInv35 := oneMinusXSqInv25 * oneMinusXSqInv 114 asinDer1DivF1Tab[i] = oneMinusXSqInv05 115 asinDer2DivF2Tab[i] = (x * oneMinusXSqInv15) * oneDivF2 116 asinDer3DivF3Tab[i] = ((1 + 2*x*x) * oneMinusXSqInv25) * oneDivF3 117 asinDer4DivF4Tab[i] = ((5 + 2*x*(2+x*(5-2*x))) * oneMinusXSqInv35) * oneDivF4 118 } 119 120 // earth radius 121 a := 6378137.0 122 b := 6356752.31420 123 a2 := a * a 124 b2 := b * b 125 earthDiameterPerLatitude = make([]float64, radiusTabsSize) 126 earthDiameterPerLatitude[0] = 2.0 * a / 1000 127 earthDiameterPerLatitude[radiusTabsSize-1] = 2.0 * b / 1000 128 for i := 1; i < radiusTabsSize-1; i++ { 129 lat := math.Pi * float64(i) / (2*radiusTabsSize - 1) 130 one := math.Pow(a2*math.Cos(lat), 2) 131 two := math.Pow(b2*math.Sin(lat), 2) 132 three := math.Pow(float64(a)*math.Cos(lat), 2) 133 four := math.Pow(b*math.Sin(lat), 2) 134 radius := math.Sqrt((one + two) / (three + four)) 135 earthDiameterPerLatitude[i] = 2 * radius / 1000 136 } 137} 138 139// earthDiameter returns an estimation of the earth's diameter at the specified 140// latitude in kilometers 141func earthDiameter(lat float64) float64 { 142 index := math.Mod(math.Abs(lat)*radiusIndexer+0.5, float64(len(earthDiameterPerLatitude))) 143 if math.IsNaN(index) { 144 return 0 145 } 146 return earthDiameterPerLatitude[int(index)] 147} 148 149var pio2 = math.Pi / 2 150 151func sin(a float64) float64 { 152 return cos(a - pio2) 153} 154 155// cos is a sloppy math (faster) implementation of math.Cos 156func cos(a float64) float64 { 157 if a < 0.0 { 158 a = -a 159 } 160 if a > sinCosMaxValueForIntModulo { 161 return math.Cos(a) 162 } 163 // index: possibly outside tables range. 164 index := int(a*sinCosIndexer + 0.5) 165 delta := (a - float64(index)*sinCosDeltaHi) - float64(index)*sinCosDeltaLo 166 // Making sure index is within tables range. 167 // Last value of each table is the same than first, so we ignore it (tabs size minus one) for modulo. 168 index &= (sinCosTabsSize - 2) // index % (SIN_COS_TABS_SIZE-1) 169 indexCos := cosTab[index] 170 indexSin := sinTab[index] 171 return indexCos + delta*(-indexSin+delta*(-indexCos*oneDivF2+delta*(indexSin*oneDivF3+delta*indexCos*oneDivF4))) 172} 173 174// asin is a sloppy math (faster) implementation of math.Asin 175func asin(a float64) float64 { 176 var negateResult bool 177 if a < 0 { 178 a = -a 179 negateResult = true 180 } 181 if a <= asinMaxValueForTabs { 182 index := int(a*asinIndexer + 0.5) 183 delta := a - float64(index)*asinDelta 184 result := asinTab[index] + delta*(asinDer1DivF1Tab[index]+delta*(asinDer2DivF2Tab[index]+delta*(asinDer3DivF3Tab[index]+delta*asinDer4DivF4Tab[index]))) 185 if negateResult { 186 return -result 187 } 188 return result 189 } 190 // value > ASIN_MAX_VALUE_FOR_TABS, or value is NaN 191 // This part is derived from fdlibm. 192 if a < 1 { 193 t := (1.0 - a) * 0.5 194 p := t * (asinPs0 + t*(asinPs1+t*(asinPs2+t*(asinPs3+t*(asinPs4+t+asinPs5))))) 195 q := 1.0 + t*(asinQs1+t*(asinQs2+t*(asinQs3+t*asinQs4))) 196 s := math.Sqrt(t) 197 z := s + s*(p/q) 198 result := asinPio2Hi - ((z + z) - asinPio2Lo) 199 if negateResult { 200 return -result 201 } 202 return result 203 } 204 // value >= 1.0, or value is NaN 205 if a == 1.0 { 206 if negateResult { 207 return -math.Pi / 2 208 } 209 return math.Pi / 2 210 } 211 return math.NaN() 212} 213