Source Code for Module gavo.utils.mathtricks

  1  """ 
  2  Math-related helper functions. 
  3  """ 
  4   
  5  #c Copyright 2008-2019, the GAVO project 
  6  #c 
  7  #c This program is free software, covered by the GNU GPL.  See the 
  8  #c COPYING file in the source distribution. 
  9   
 10   
 11  import math 
 12   
 13  from gavo.utils import codetricks 
 14   
 15  DEG = math.pi/180 
 16  ARCSEC = DEG/3600 
 17   
 18 -def findMinimum(f, left, right, minInterval=3e-8): 
 19          """returns an estimate for the minimum of the single-argument function f  
 20          on (left,right). 
 21   
 22          minInterval is a fourth of the smallest test interval considered.   
 23   
 24          For constant functions, a value close to left will be returned. 
 25   
 26          This function should only be used on functions having exactly 
 27          one minimum in the interval. 
 28          """ 
 29  # replace this at some point by some better method (Num. Recip. in C, 394f) 
 30  # -- this is easy to fool and massively suboptimal. 
 31          mid = (right+left)/2. 
 32          offset = (right-left)/4. 
 33          if offset<minInterval: 
 34                  return mid 
 35          if f(left+offset)<=f(mid+offset): 
 36                  return findMinimum(f, left, mid, minInterval) 
 37          else: 
 38                  return findMinimum(f, mid, right, minInterval) 
 39   
 40   
 41 -class getHexToBin(codetricks.CachedResource): 
 42          """returns a dictionary mapping hex chars to their binary expansions. 
 43          """ 
 44          @classmethod 
 45 -        def impl(cls): 
 46                  return dict(zip( 
 47                          "0123456789abcdef", 
 48                          ["0000", "0001", "0010", "0011", "0100", "0101", "0110", "0111", 
 49                           "1000", "1001", "1010", "1011", "1100", "1101", "1110", "1111",])) 
 50                   
 51   
 52 -def toBinary(anInt, desiredLength=None): 
 53          """returns anInt as a string with its binary digits, MSB first. 
 54   
 55          If desiredLength is given and the binary expansion is shorter, 
 56          the value will be padded with zeros. 
 57   
 58          >>> toBinary(349) 
 59          '101011101' 
 60          >>> toBinary(349, 10) 
 61          '0101011101' 
 62          """ 
 63          h2b = getHexToBin() 
 64          res = "".join(h2b[c] for c in "%x"%anInt).lstrip("0") 
 65          if desiredLength is not None: 
 66                  res = "0"*(desiredLength-len(res))+res 
 67          return res 
 68   
 69   
 70 -def _dotprod3(seq1, seq2, inds=(0,1,2)): 
 71          return sum(seq1[i]*seq2[i] for i in inds) 
 72   
 73   
 74 -class Matrix3(object): 
 75          """A quick and easy 3d matrix. 
 76   
 77          This is just so we don't depend on numpy for trivial stuff.  The 
 78          components are stored in a tuple of rows. 
 79          """ 
 80          indices = range(3) 
 81   
 82 -        def __init__(self, row1, row2, row3): 
 83                  self.rows = (tuple(row1), tuple(row2), tuple(row3)) 
 84   
 85 -        def __eq__(self, other): 
 86                  return (isinstance(other, Matrix3) 
 87                          and self.rows==other.rows) 
 88           
 89 -        def __ne__(self, other): 
 90                  return not self.__eq__(other) 
 91   
 92 -        def vecMul(self, vec): 
 93                  """returns the result of right-multiplying self to vec. 
 94   
 95                  The sequence vec is interpreted as a column vector. 
 96                  """ 
 97                  return tuple(_dotprod3(self.rows[i], vec) for i in self.indices) 
 98   
 99 -        def matMul(self, mat): 
100                  """returns the result of multiplying mat to self from the right. 
101                  """ 
102                  cols = mat.getColumns() 
103                  return self.__class__(*tuple( 
104                                  tuple(_dotprod3(row, col) for col in cols) 
105                          for row in self.rows)) 
106           
107 -        def getColumns(self): 
108                  """returns the column vectors of this matrix in a 3-tuple. 
109                  """ 
110                  return tuple( 
111                                  tuple(self.rows[rowInd][colInd] for rowInd in self.indices) 
112                          for colInd in self.indices) 
113   
114   
115 -def _normalize3(vec): 
116          len = math.sqrt(_dotprod3(vec, vec)) 
117          return tuple(c/len for c in vec) 
118   
119   
120 -def getRotX(angle): 
121          """returns a 3-rotation matrix for rotating angle radians around x. 
122          """ 
123          c, s = math.cos(angle), math.sin(angle) 
124          return Matrix3((1, 0, 0), (0, c, s), (0, -s, c)) 
125   
126   
127 -def getRotZ(angle): 
128          """returns a 3-rotation matrix for rotating angle radians around z. 
129          """ 
130          c, s = math.cos(angle), math.sin(angle) 
131          return Matrix3((c, s, 0), (-s, c, 0), (0, 0, 1)) 
132   
133   
134 -def spherToCart(theta, phi): 
135          """returns a 3-cartesian unit vector pointing to longitude theta, 
136          latitude phi. 
137   
138          The angles are in rad. 
139          """ 
140          cp = math.cos(phi) 
141          return math.cos(theta)*cp, math.sin(theta)*cp, math.sin(phi) 
142   
143   
144 -def cartToSpher(unitvector): 
145          """returns spherical coordinates for a 3-unit vector. 
146   
147          We do not check if unitvector actually *is* a unit vector.  The returned 
148          angles are in rad. 
149          """ 
150          x, y, z = unitvector 
151          rInXY = math.sqrt(x**2+y**2) 
152          if abs(rInXY)<1e-9:  # pole 
153                  theta = 0 
154          else: 
155                  theta = math.atan2(y, x) 
156          if theta<0: 
157                  theta += 2*math.pi 
158          phi = math.atan2(z, rInXY) 
159          return (theta, phi) 
160   
161   
162 -def spherDist(vec1, vec2): 
163          """returns the spherical distance (in radian) between the unit vectors 
164          vec1 and vec2. 
165          """ 
166          return math.acos(_dotprod3(vec1, vec2)) 
167   
168   
169 -def _test(): 
170          import doctest, mathtricks 
171          doctest.testmod(mathtricks) 
172   
173   
174  if __name__=="__main__": 
175          _test() 
176