Source code for gavo.utils.mathtricks

"""
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()