213 lines
6.2 KiB
Python
213 lines
6.2 KiB
Python
"""geomath.py: transcription of GeographicLib::Math class."""
|
|
# geomath.py
|
|
#
|
|
# This is a rather literal translation of the GeographicLib::Math class to
|
|
# python. See the documentation for the C++ class for more information at
|
|
#
|
|
# https://geographiclib.sourceforge.io/html/annotated.html
|
|
#
|
|
# Copyright (c) Charles Karney (2011-2019) <charles@karney.com> and
|
|
# licensed under the MIT/X11 License. For more information, see
|
|
# https://geographiclib.sourceforge.io/
|
|
######################################################################
|
|
|
|
import sys
|
|
import math
|
|
|
|
class Math(object):
|
|
"""
|
|
Additional math routines for GeographicLib.
|
|
|
|
This defines constants:
|
|
epsilon, difference between 1 and the next bigger number
|
|
digits, the number of digits in the fraction of a real number
|
|
minval, minimum normalized positive number
|
|
maxval, maximum finite number
|
|
nan, not a number
|
|
inf, infinity
|
|
"""
|
|
|
|
digits = 53
|
|
epsilon = math.pow(2.0, 1-digits)
|
|
minval = math.pow(2.0, -1022)
|
|
maxval = math.pow(2.0, 1023) * (2 - epsilon)
|
|
inf = float("inf") if sys.version_info > (2, 6) else 2 * maxval
|
|
nan = float("nan") if sys.version_info > (2, 6) else inf - inf
|
|
|
|
def sq(x):
|
|
"""Square a number"""
|
|
|
|
return x * x
|
|
sq = staticmethod(sq)
|
|
|
|
def cbrt(x):
|
|
"""Real cube root of a number"""
|
|
|
|
y = math.pow(abs(x), 1/3.0)
|
|
return y if x > 0 else (-y if x < 0 else x)
|
|
cbrt = staticmethod(cbrt)
|
|
|
|
def log1p(x):
|
|
"""log(1 + x) accurate for small x (missing from python 2.5.2)"""
|
|
|
|
if sys.version_info > (2, 6):
|
|
return math.log1p(x)
|
|
|
|
y = 1 + x
|
|
z = y - 1
|
|
# Here's the explanation for this magic: y = 1 + z, exactly, and z
|
|
# approx x, thus log(y)/z (which is nearly constant near z = 0) returns
|
|
# a good approximation to the true log(1 + x)/x. The multiplication x *
|
|
# (log(y)/z) introduces little additional error.
|
|
return x if z == 0 else x * math.log(y) / z
|
|
log1p = staticmethod(log1p)
|
|
|
|
def atanh(x):
|
|
"""atanh(x) (missing from python 2.5.2)"""
|
|
|
|
if sys.version_info > (2, 6):
|
|
return math.atanh(x)
|
|
|
|
y = abs(x) # Enforce odd parity
|
|
y = Math.log1p(2 * y/(1 - y))/2
|
|
return y if x > 0 else (-y if x < 0 else x)
|
|
atanh = staticmethod(atanh)
|
|
|
|
def copysign(x, y):
|
|
"""return x with the sign of y (missing from python 2.5.2)"""
|
|
|
|
if sys.version_info > (2, 6):
|
|
return math.copysign(x, y)
|
|
|
|
return math.fabs(x) * (-1 if y < 0 or (y == 0 and 1/y < 0) else 1)
|
|
copysign = staticmethod(copysign)
|
|
|
|
def norm(x, y):
|
|
"""Private: Normalize a two-vector."""
|
|
r = math.hypot(x, y)
|
|
return x/r, y/r
|
|
norm = staticmethod(norm)
|
|
|
|
def sum(u, v):
|
|
"""Error free transformation of a sum."""
|
|
# Error free transformation of a sum. Note that t can be the same as one
|
|
# of the first two arguments.
|
|
s = u + v
|
|
up = s - v
|
|
vpp = s - up
|
|
up -= u
|
|
vpp -= v
|
|
t = -(up + vpp)
|
|
# u + v = s + t
|
|
# = round(u + v) + t
|
|
return s, t
|
|
sum = staticmethod(sum)
|
|
|
|
def polyval(N, p, s, x):
|
|
"""Evaluate a polynomial."""
|
|
y = float(0 if N < 0 else p[s]) # make sure the returned value is a float
|
|
while N > 0:
|
|
N -= 1; s += 1
|
|
y = y * x + p[s]
|
|
return y
|
|
polyval = staticmethod(polyval)
|
|
|
|
def AngRound(x):
|
|
"""Private: Round an angle so that small values underflow to zero."""
|
|
# The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
|
|
# for reals = 0.7 pm on the earth if x is an angle in degrees. (This
|
|
# is about 1000 times more resolution than we get with angles around 90
|
|
# degrees.) We use this to avoid having to deal with near singular
|
|
# cases when x is non-zero but tiny (e.g., 1.0e-200).
|
|
z = 1/16.0
|
|
y = abs(x)
|
|
# The compiler mustn't "simplify" z - (z - y) to y
|
|
if y < z: y = z - (z - y)
|
|
return 0.0 if x == 0 else (-y if x < 0 else y)
|
|
AngRound = staticmethod(AngRound)
|
|
|
|
def remainder(x, y):
|
|
"""remainder of x/y in the range [-y/2, y/2]."""
|
|
z = math.fmod(x, y) if Math.isfinite(x) else Math.nan
|
|
# On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) = +0.0
|
|
# This fixes this bug. See also Math::AngNormalize in the C++ library.
|
|
# sincosd has a similar fix.
|
|
z = x if x == 0 else z
|
|
return (z + y if z < -y/2 else
|
|
(z if z < y/2 else z -y))
|
|
remainder = staticmethod(remainder)
|
|
|
|
def AngNormalize(x):
|
|
"""reduce angle to (-180,180]"""
|
|
|
|
y = Math.remainder(x, 360)
|
|
return 180 if y == -180 else y
|
|
AngNormalize = staticmethod(AngNormalize)
|
|
|
|
def LatFix(x):
|
|
"""replace angles outside [-90,90] by NaN"""
|
|
|
|
return Math.nan if abs(x) > 90 else x
|
|
LatFix = staticmethod(LatFix)
|
|
|
|
def AngDiff(x, y):
|
|
"""compute y - x and reduce to [-180,180] accurately"""
|
|
|
|
d, t = Math.sum(Math.AngNormalize(-x), Math.AngNormalize(y))
|
|
d = Math.AngNormalize(d)
|
|
return Math.sum(-180 if d == 180 and t > 0 else d, t)
|
|
AngDiff = staticmethod(AngDiff)
|
|
|
|
def sincosd(x):
|
|
"""Compute sine and cosine of x in degrees."""
|
|
|
|
r = math.fmod(x, 360) if Math.isfinite(x) else Math.nan
|
|
q = 0 if Math.isnan(r) else int(round(r / 90))
|
|
r -= 90 * q; r = math.radians(r)
|
|
s = math.sin(r); c = math.cos(r)
|
|
q = q % 4
|
|
if q == 1:
|
|
s, c = c, -s
|
|
elif q == 2:
|
|
s, c = -s, -c
|
|
elif q == 3:
|
|
s, c = -c, s
|
|
# Remove the minus sign on -0.0 except for sin(-0.0).
|
|
# On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) = +0.0
|
|
# (x, c) here fixes this bug. See also Math::sincosd in the C++ library.
|
|
# AngNormalize has a similar fix.
|
|
s, c = (x, c) if x == 0 else (0.0+s, 0.0+c)
|
|
return s, c
|
|
sincosd = staticmethod(sincosd)
|
|
|
|
def atan2d(y, x):
|
|
"""compute atan2(y, x) with the result in degrees"""
|
|
|
|
if abs(y) > abs(x):
|
|
q = 2; x, y = y, x
|
|
else:
|
|
q = 0
|
|
if x < 0:
|
|
q += 1; x = -x
|
|
ang = math.degrees(math.atan2(y, x))
|
|
if q == 1:
|
|
ang = (180 if y >= 0 else -180) - ang
|
|
elif q == 2:
|
|
ang = 90 - ang
|
|
elif q == 3:
|
|
ang = -90 + ang
|
|
return ang
|
|
atan2d = staticmethod(atan2d)
|
|
|
|
def isfinite(x):
|
|
"""Test for finiteness"""
|
|
|
|
return abs(x) <= Math.maxval
|
|
isfinite = staticmethod(isfinite)
|
|
|
|
def isnan(x):
|
|
"""Test if nan"""
|
|
|
|
return math.isnan(x) if sys.version_info > (2, 6) else x != x
|
|
isnan = staticmethod(isnan)
|