"""
Math-related helper functions.
"""
#c Copyright 2008-2023, the GAVO project <gavo@ari.uni-heidelberg.de>
#c
#c This program is free software, covered by the GNU GPL. See the
#c COPYING file in the source distribution.
import functools
import math
from astropy.coordinates import ( # type: ignore
UnitSphericalRepresentation, CartesianRepresentation)
from astropy import units as u # type: ignore
import numpy
from gavo.utils.dachstypes import (Callable, Dict, Optional)
DEG = math.pi/180
ARCSEC = DEG/3600
[docs]def findMinimum(
f: Callable[[float], float],
left: float,
right: float,
minInterval: float=3e-8) -> float:
"""returns an estimate for the minimum of the single-argument function f
on (left,right).
minInterval is a fourth of the smallest test interval considered.
For constant functions, a value close to left will be returned.
This function should only be used on functions having exactly
one minimum in the interval.
"""
# replace this at some point by some better method (Num. Recip. in C, 394f)
# -- this is easy to fool and massively suboptimal.
mid = (right+left)/2.
offset = (right-left)/4.
if offset<minInterval:
return mid
if f(left+offset)<=f(mid+offset):
return findMinimum(f, left, mid, minInterval)
else:
return findMinimum(f, mid, right, minInterval)
[docs]@functools.cache
def getHexToBin() -> Dict[str, str]:
"""returns a dictionary mapping hex chars to their binary expansions.
"""
return dict(list(zip(
"0123456789abcdef",
["0000", "0001", "0010", "0011", "0100", "0101", "0110", "0111",
"1000", "1001", "1010", "1011", "1100", "1101", "1110", "1111",])))
[docs]def toBinary(anInt: int, desiredLength: Optional[int]=None) -> str:
"""returns anInt as a string with its binary digits, MSB first.
If desiredLength is given and the binary expansion is shorter,
the value will be padded with zeros.
>>> toBinary(349)
'101011101'
>>> toBinary(349, 10)
'0101011101'
"""
h2b = getHexToBin()
res = "".join(h2b[c] for c in "%x"%anInt).lstrip("0")
if desiredLength is not None:
res = "0"*(desiredLength-len(res))+res
return res
[docs]def roundO2M(num: float) -> int:
"""returns a plausible rounding of num.
This will round up the last couple of digits. For now, this will only
do cardinals.
>>> roundO2M(0)
0
>>> roundO2M(2.5)
2
>>> roundO2M(15)
20
>>> roundO2M(9900)
10000
>>> roundO2M(8321)
8400
>>> roundO2M(3.2349302e9)
3300000000
"""
if num<10:
return int(round(num))
mant = ("%f"%num).split(".")[0]
if len(mant)==2:
nextVal = int(mant[0])+1
return int("%s0"%nextVal)
nextVal = int(mant[:2])+1
return int("%s%s"%(nextVal, "0"*(len(mant)-2)))
[docs]def getRotX(angle):
"""returns a 3-rotation matrix for rotating angle radians around x.
"""
c, s = math.cos(angle), math.sin(angle)
return numpy.array([(1, 0, 0), (0, c, s), (0, -s, c)])
[docs]def getRotY(angle):
"""returns a 3-rotation matrix for rotating angle radians around y.
"""
c, s = math.cos(angle), math.sin(angle)
return numpy.array([[c, 0, -s], [0, 1, 0], [s, 0, c]])
[docs]def getRotZ(angle):
"""returns a 3-rotation matrix for rotating angle radians around z.
"""
c, s = math.cos(angle), math.sin(angle)
return numpy.array([(c, s, 0), (-s, c, 0), (0, 0, 1)])
[docs]def spherToCart(theta, phi):
"""returns a 3-cartesian unit vector pointing to longitude theta,
latitude phi.
The angles are in rad.
"""
return UnitSphericalRepresentation(theta*u.rad, phi*u.rad
).represent_as(CartesianRepresentation).xyz.value
[docs]def cartToSpher(unitvector):
"""returns spherical coordinates for a 3-unit vector.
We do not check if unitvector actually *is* a unit vector. The returned
angles are in rad.
"""
spherical = CartesianRepresentation(*unitvector
).represent_as(UnitSphericalRepresentation)
return [spherical.lon.to(u.rad).value, spherical.lat.to(u.rad).value]
[docs]def floorInt(flt: float) -> float:
"""returns the largest integer smaller than flt, except if flt is a
special float, in which case flt is returned.
"""
if math.isfinite(flt):
return int(math.floor(flt))
return flt
[docs]def ceilInt(flt: float) -> float:
"""returns the smallest integer largest than flt, except if flt is a
special float, in which case flt is returned.
"""
if math.isfinite(flt):
return int(math.ceil(flt))
return flt
if __name__=="__main__": # pragma: no cover
import doctest
doctest.testmod()