Source Code for Module gavo.stc.spherc

  1  """ 
  2  Spherical sky coordinates and helpers. 
  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  # XXX TODO: Replace the actual transformations performed here with 
 12  # stuff from astropy. 
 13   
 14  from math import sin, cos 
 15  import math 
 16   
 17  import numpy 
 18  from numpy import linalg as la 
 19   
 20  from gavo.stc import common 
 21  from gavo.stc import sphermath 
 22  from gavo.stc import times 
 23  from gavo.stc import units 
 24  from gavo.utils import DEG, ARCSEC, memoized 
 25   
 26   
 27  ############### Basic definitions for transforms 
 28   
 29  # Finding transformation sequences: This, in principle, is a standard 
 30  # graph problem.  However, we have lots of underspecified transforms, 
 31  # which makes building a Dijkstrable graph somewhat inconvenient.  So, 
 32  # instead of striving for an optimal shortest path, we go for a 
 33  # greedy search with some heuristics.  The most nonstandard feature is 
 34  # that nodes are built ad-hoc and noncircularity of the paths thorugh 
 35  # the "virtual" graph is checked on these ad-hoc nodes. 
 36   
 37  # The nodes in the graph are triples of (system, equinox, refpos).  The 
 38  # vertices are triples (fromNode, toNode, transform generator). 
 39  # Transform generators are functions  
 40  # 
 41  # f(fromNode, toNode) -> (function or matrix) 
 42  # 
 43  # The arguments are node triples, the result either a function taking 
 44  # and returning 6-vectors or numpy matrices.  These functions may 
 45  # assume that only "appropriate" values are passed in as nodes, i.e., 
 46  # they are not assumed to check that the are actually able to produce 
 47  # the requested transformation. 
 48   
 49 -class _Wildcard(object):  
 50          """is an object that compares equal to everything. 
 51   
 52          This is used for underspecification of transforms, see SAME and 
 53          ANYVAL below. 
 54          """ 
 55 -        def __init__(self, name): 
 56                  self.name = name 
 57           
 58 -        def __repr__(self): 
 59                  return self.name 
 60   
 61 -        def __ne__(self, other): return False 
 62 -        def __eq__(self, other): return True 
 63   
 64   
 65  SAME = _Wildcard("SAME") 
 66  ANYVAL = _Wildcard("ANYVAL") 
 67   
 68 -def _specifyTrafo(trafo, fromTuple, toTuple): 
 69          """fills in underspecified values in trafo from fromTuple and toTuple. 
 70   
 71          The rules are: In the source, both SAME and ANYVAL are filled from 
 72          fromTuple.  In the destination, SAME is filled from fromTuple, 
 73          ANYVAL is filled from toTuple. 
 74   
 75          The function returns a new transformation triple. 
 76          """ 
 77          src, dst, tgen = trafo 
 78          newSrc, newDst = [], [] 
 79          for ind, val in enumerate(src): 
 80                  if val is SAME or val is ANYVAL: 
 81                          newSrc.append(fromTuple[ind]) 
 82                  else: 
 83                          newSrc.append(val) 
 84          for ind, val in enumerate(dst): 
 85                  if val is SAME: 
 86                          newDst.append(fromTuple[ind]) 
 87                  elif val is ANYVAL: 
 88                          newDst.append(toTuple[ind]) 
 89                  else: 
 90                          newDst.append(val) 
 91          return tuple(newSrc), tuple(newDst), tgen 
 92   
 93   
 94 -def _makeFindPath(transforms): 
 95          """returns a function for path finding in the virtual graph 
 96          defined by transforms. 
 97   
 98          Each transform is a triple of (fromNode, toNode, transformFactory). 
 99   
100          There's quite a bit of application-specific heuristics built in 
101          here, so there's litte you can do with this code outside of 
102          STC transforms construction. 
103          """ 
104          def findPath(fromTuple, toTuple, path=()): 
105                  """returns a sequence of transformation triples that lead from 
106                  fromTuple to toTuple. 
107   
108                  fromTuple and toTuple are node triples (i.e., (system, equinox, 
109                  refpoint)). 
110   
111                  The returned path is not guaranteed to be the shortest or even 
112                  the numerically most stable.  It is the result of a greedy 
113                  search for a cycle free path between the two "non-virtual" nodes. 
114                  To keep the paths reasonable, we apply the heuristic that 
115                  transformations keeping the system are preferable. 
116   
117                  The simple heuristics sometimes need help; e.g., the transformations 
118                  below add explicit transformations to j2000 and b1950; you will always 
119                  need this if your transformations include "magic" values for otherwise 
120                  underspecified items. 
121                  """ 
122                  seenSystems = set(c[0] for c in path) | set(c[1] for c in path) 
123                  candidates = [_specifyTrafo(t, fromTuple, toTuple)  
124                                  for t in transforms if t[0]==fromTuple] 
125                  # sort operations within the same reference system to the start 
126                  candidates = [c for c in candidates if c[1][0]==toTuple[0]] + [ 
127                          c for c in candidates if c[1]!=toTuple[0]] 
128                  for cand in candidates: 
129                          srcSystem, dstSystem, tgen = cand 
130                          # Ignore identities or trafos leading to cycles 
131                          if srcSystem==dstSystem or dstSystem in seenSystems: 
132                                  continue 
133                          if dstSystem==toTuple:  # If we are done, return result 
134                                  return path+(cand,) 
135                          else: 
136                                  # Do the depth-first search 
137                                  np = findPath(dstSystem, toTuple, path+(cand,)) 
138                                  if np: 
139                                          return np 
140          return findPath 
141   
142   
143 -def tupleMin(t1, t2): 
144          """returns the element-wise minimum of two tuples: 
145          """ 
146          return tuple(min(i1, i2) for i1, i2 in zip(t1, t2)) 
147   
148   
149 -def tupleMax(t1, t2): 
150          """returns the element-wise maximum of two tuples: 
151          """ 
152          return tuple(max(i1, i2) for i1, i2 in zip(t1, t2)) 
153   
154   
155  ############### Computation of precession matrices. 
156   
157 -def prec_IAU1976(fromEquinox, toEquinox): 
158          """returns the precession angles in the IAU 1976 system. 
159   
160          The expressions are those of Lieske, A&A 73, 282. 
161   
162          This function is for the precTheory argument of getPrecMatrix. 
163          """ 
164          # time differences have to be in julian centuries 
165          # captial T in Lieske 
166          fromDiff = times.getSeconds(fromEquinox-times.dtJ2000)/common.secsPerJCy  
167          # lowercase T in Lieske 
168          toDiff = times.getSeconds(toEquinox-fromEquinox)/common.secsPerJCy   
169   
170          # Lieske's expressions yield arcsecs, fix below 
171          zeta = toDiff*(2306.2181+1.39656*fromDiff-0.000139*fromDiff**2 
172                  )+toDiff**2*(0.30188-0.000344*fromDiff 
173                  )+toDiff**3*0.017998 
174          z =    toDiff*(2306.2181+1.39656*fromDiff-0.000139*fromDiff**2 
175                  )+toDiff**2*(1.09468+0.000066*fromDiff 
176                  )+toDiff**3*0.018203 
177          theta = toDiff*(2004.3109-0.85330*fromDiff-0.000217*fromDiff**2 
178                  )-toDiff**2*(0.42665+0.000217*fromDiff 
179                  )-toDiff**3*0.041833 
180          return zeta*ARCSEC, z*ARCSEC, theta*ARCSEC 
181   
182   
183  _dtB1850 = times.bYearToDateTime(1850) 
184   
185 -def prec_Newcomb(fromEquinox, toEquinox): 
186          """returns the precession angles for the newcomp 
187   
188          This function is for the precTheory argument of getPrecMatrix. 
189   
190          The expressions are Kinoshita's (1975)'s (SAOSR 364)  
191          This is somewhat at odds with Yallop's choice of Andoyer in the FK4-FK5 
192          machinery below, but that really shouldn't matter. 
193          """ 
194          # time differences have to be in tropical centuries 
195          T = times.getSeconds(fromEquinox-_dtB1850)/(common.tropicalYear*86400*100) 
196          t = times.getSeconds(toEquinox-fromEquinox)/(common.tropicalYear*86400*100) 
197   
198          polyVal = 2303.5548+(1.39720+0.000059*T)*T 
199          zeta = (polyVal+(0.30242-0.000269*T+0.017996*t)*t)*t 
200          z = (polyVal+(1.09478+0.000387*T+0.018324*t)*t)*t 
201          theta = (2005.1125+(-0.85294-0.000365*T)*T 
202                  +(-0.42647-0.000365*T-0.041802*t)*t)*t 
203          return zeta*ARCSEC, z*ARCSEC, theta*ARCSEC 
204   
205   
206 -def getPrecMatrix(fromEquinox, toEquinox, precTheory): 
207          """returns a precession matrix in the IAU(1976) system. 
208   
209          fromEquinox and toEquinox are datetimes (in case of doubt, TT). 
210   
211          precTheory(fromEquinox, toEquinox) -> zeta, z, theta computes the 
212          precession angles.  Defined in this module are prec_IAU1976 and  
213          prec_Newcomb, but you can provide your own.  The angles must all be 
214          in rad. 
215          """ 
216          zeta, z, theta = precTheory(fromEquinox, toEquinox) 
217          return numpy.dot( 
218                  numpy.dot(sphermath.getRotZ(-z), sphermath.getRotY(theta)), 
219                  sphermath.getRotZ(-zeta)) 
220   
221   
222  _nullMatrix = numpy.zeros((3,3)) 
223 -def threeToSix(matrix): 
224          """returns a 6-matrix from a 3-matrix suitable for precessing our 
225          6-vectors. 
226          """ 
227          return numpy.concatenate(( 
228                  numpy.concatenate(  (matrix,      _nullMatrix), 1), 
229                  numpy.concatenate(  (_nullMatrix, matrix     ), 1))) 
230   
231 -def _getFullPrecMatrix(fromNode, toNode, precTheory): 
232          """returns a full 6x6 matrix for transforming positions and proper motions. 
233   
234          This only works for proper equatorial coordinates in both STC values. 
235   
236          precTheory is a function returning precession angles. 
237          """ 
238          return threeToSix(getPrecMatrix(fromNode[1], toNode[1], precTheory)) 
239   
240   
241 -def _getNewcombPrecMatrix(fromNode, toNode, sixTrans): 
242          return _getFullPrecMatrix(fromNode, toNode, prec_Newcomb) 
243   
244 -def _getIAU1976PrecMatrix(fromNode, toNode, sixTrans): 
245          return _getFullPrecMatrix(fromNode, toNode, prec_IAU1976) 
246   
247   
248  ############### FK4-FK5 system transformation 
249  # This follows the prescription of Yallop et al, AJ 97, 274 
250   
251  # Transformation matrix according to Yallop 
252  _fk4ToFK5MatrixYallop = numpy.array([ 
253          [0.999925678186902, -0.011182059642247, -0.004857946558960, 
254                  0.000002423950176, -0.000000027106627, -0.000000011776558], 
255          [0.011182059571766, 0.999937478448132, -0.000027176441185, 
256                  0.000000027106627, 0.000002723978783, -0.000000000065874], 
257          [0.004857946721186, -0.000027147426489, 0.9999881997387700, 
258                  0.000000011776559, -0.000000000065816, 0.000002424101735], 
259          [-0.000541652366951, -0.237968129744288, 0.436227555856097, 
260                  0.999947035154614, -0.011182506121805, -0.004857669684959], 
261          [0.237917612131583, -0.002660763319071, -0.008537771074048, 
262                  0.011182506007242, 0.999958833818833, -0.000027184471371], 
263          [-0.436111276039270, 0.012259092261564, 0.002119110818172, 
264                  0.004857669948650, -0.000027137309539, 1.000009560363559]]) 
265   
266  # Transformation matrix according to SLALIB-F 
267  _fk4ToFK5MatrixSla = numpy.transpose(numpy.array([ 
268          [+0.9999256782, +0.0111820610, +0.0048579479,  
269                  -0.000551, +0.238514, -0.435623], 
270          [-0.0111820611, +0.9999374784, -0.0000271474,  
271                  -0.238565, -0.002667, +0.012254], 
272          [-0.0048579477, -0.0000271765, +0.9999881997,  
273                  +0.435739, -0.008541, +0.002117], 
274          [+0.00000242395018, +0.00000002710663, +0.00000001177656,  
275                  +0.99994704, +0.01118251, +0.00485767], 
276          [-0.00000002710663, +0.00000242397878, -0.00000000006582,  
277                  -0.01118251, +0.99995883, -0.00002714], 
278          [-0.00000001177656, -0.00000000006587, +0.00000242410173,  
279                  -0.00485767, -0.00002718, +1.00000956]])) 
280   
281  # Inverse transformation matrix according to SLALIB-F 
282  _fk5ToFK4Matrix = numpy.transpose(numpy.array([ 
283  [+0.9999256795, -0.0111814828, -0.0048590040,  
284          -0.000551, -0.238560, +0.435730], 
285  [+0.0111814828, +0.9999374849, -0.0000271557,  
286          +0.238509, -0.002667, -0.008541], 
287  [+0.0048590039, -0.0000271771, +0.9999881946,  
288          -0.435614, +0.012254, +0.002117], 
289  [-0.00000242389840, +0.00000002710544, +0.00000001177742,  
290          +0.99990432, -0.01118145, -0.00485852], 
291  [-0.00000002710544, -0.00000242392702, +0.00000000006585,  
292          +0.01118145, +0.99991613, -0.00002716], 
293  [-0.00000001177742, +0.00000000006585, -0.00000242404995,  
294          +0.00485852, -0.00002717, +0.99996684]])) 
295   
296   
297   
298  # Positional correction due to E-Terms, in rad (per tropical century in the 
299  # case of Adot, which is ok for sphermath._svPosUnit (Yallop et al, loc cit, p. 
300  # 276)).  We ignore the difference between tropical and julian centuries. 
301  _b1950ETermsPos = numpy.array([-1.62557e-6, -0.31919e-6, -0.13843e-6]) 
302  _b1950ETermsVel = numpy.array([1.245e-3, -1.580e-3, -0.659e-3]) 
303  _yallopK = common.secsPerJCy/(units.oneAU/1e3) 
304  _yallopKSla = 21.095; 
305  _pcPerCyToKmPerSec = units.getRedshiftConverter("pc", "cy", "km", "s") 
306   
307  # An SVConverter to bring 6-vectors to the spherical units Yallop prescribes. 
308  _yallopSVConverter = sphermath.SVConverter((0,0,0), ("rad", "rad", "arcsec"),  
309          (0,0,0), ("arcsec", "arcsec", "km"), ("cy", "cy", "s")) 
310   
311   
312 -def _svToYallop(sv, yallopK): 
313          """returns r and rdot vectors suitable for Yallop's recipe from our 
314          6-vectors. 
315          """ 
316          (alpha, delta, prlx), (pma, pmd, rv) = _yallopSVConverter.from6(sv) 
317   
318          yallopR = numpy.array([cos(alpha)*cos(delta), 
319                  sin(alpha)*cos(delta), sin(delta)]) 
320          yallopRd = numpy.array([ 
321                  -pma*sin(alpha)*cos(delta)-pmd*cos(alpha)*sin(delta), 
322                  pma*cos(alpha)*cos(delta)-pmd*sin(alpha)*sin(delta), 
323                  pmd*cos(delta)])+yallopK*rv*prlx*yallopR 
324          return yallopR, yallopRd, (rv, prlx) 
325   
326   
327 -def _yallopToSv(yallop6, yallopK, rvAndPrlx): 
328          """returns a 6-Vector from a yallop-6 vector. 
329   
330          rvAndPrlx is the third item of the return value of _svToYallop. 
331          """ 
332          rv, prlx = rvAndPrlx 
333          x,y,z,xd,yd,zd = yallop6 
334          rxy2 = x**2+y**2 
335          r = math.sqrt(z**2+rxy2) 
336          if rxy2==0: 
337                  raise common.STCValueError("No spherical proper motion on poles.") 
338          alpha = math.atan2(y, x) 
339          if alpha<0: 
340                  alpha += 2*math.pi 
341          delta = math.atan2(z, math.sqrt(rxy2)) 
342          pma = (x*yd-y*xd)/rxy2 
343          pmd = (zd*rxy2-z*(x*xd+y*yd))/r/r/math.sqrt(rxy2) 
344          if abs(prlx)>1/sphermath.defaultDistance: 
345                  rv = numpy.dot(yallop6[:3], yallop6[3:])/yallopK/prlx/r 
346                  prlx = prlx/r 
347          return _yallopSVConverter.to6((alpha, delta, prlx), (pma, pmd, rv)) 
348   
349   
350 -def fk4ToFK5(sixTrans, svfk4): 
351          """returns an FK5 2000 6-vector for an FK4 1950 6-vector. 
352   
353          The procedure used is described in Yallop et al, AJ 97, 274.  E-terms 
354          of aberration are always removed from proper motions, regardless of 
355          whether the objects are within 10 deg of the pole. 
356          """ 
357          if sixTrans.slaComp: 
358                  transMatrix = _fk4ToFK5MatrixSla 
359                  yallopK = _yallopKSla 
360          else: 
361                  transMatrix = _fk4ToFK5MatrixYallop 
362                  yallopK = _yallopK 
363          yallopR, yallopRd, rvAndPrlx = _svToYallop(svfk4, yallopK) 
364   
365          # Yallop's recipe starts here 
366          if not sixTrans.slaComp:  # include Yallop's "small terms" in PM 
367                  yallopVE = (yallopRd-_b1950ETermsVel 
368                          +numpy.dot(yallopR, _b1950ETermsVel)*yallopR 
369                          +numpy.dot(yallopRd, _b1950ETermsPos)*yallopR 
370                          +numpy.dot(yallopRd, _b1950ETermsPos)*yallopRd) 
371          else: 
372                  yallopVE = (yallopRd-_b1950ETermsVel 
373                          +numpy.dot(yallopR, _b1950ETermsVel)*yallopR) 
374   
375          yallop6 = numpy.concatenate((yallopR-(_b1950ETermsPos- 
376                          numpy.dot(yallopR, _b1950ETermsPos)*yallopR), 
377                  yallopVE)) 
378          cnv = numpy.dot(transMatrix, yallop6) 
379          return _yallopToSv(cnv, yallopK, rvAndPrlx) 
380   
381   
382 -def fk5ToFK4(sixTrans, svfk5): 
383          """returns an FK4 1950 6-vector for an FK5 2000 6-vector. 
384   
385          This is basically a reversal of fk4ToFK5, except we're always operating 
386          in slaComp mode here. 
387   
388          """ 
389          yallopR, yallopRd, rvAndPrlx = _svToYallop(svfk5, _yallopKSla) 
390   
391          # first apply rotation... 
392          cnv = numpy.dot(_fk5ToFK4Matrix,  
393                  numpy.concatenate((yallopR, yallopRd))) 
394          # ... then handle E-Terms; direct inversion of Yallop's equations is 
395          # troublesome, so I basically follow what slalib does. 
396          yallopR, yallopRd = cnv[:3], cnv[3:] 
397          spatialCorr = numpy.dot(yallopR, _b1950ETermsPos)*yallopR 
398          newRMod = sphermath.vabs(yallopR+_b1950ETermsPos* 
399                  sphermath.vabs(yallopR)-spatialCorr) 
400          newR = yallopR+_b1950ETermsPos*newRMod-spatialCorr 
401          newRd = yallopRd+_b1950ETermsVel*newRMod-numpy.dot( 
402                  yallopR, _b1950ETermsVel)*yallopR 
403   
404          return _yallopToSv(numpy.concatenate((newR, newRd)), 
405                  _yallopKSla, rvAndPrlx) 
406   
407   
408  ############### Galactic coordinates 
409   
410  _galB1950pole = (192.25*DEG, 27.4*DEG) 
411  _galB1950zero = (265.6108440311*DEG, -28.9167903484*DEG) 
412   
413  _b1950ToGalTrafo = sphermath.computeTransMatrixFromPole( 
414          _galB1950pole, _galB1950zero) 
415  _b1950ToGalMatrix = threeToSix(_b1950ToGalTrafo) 
416   
417  # For convenience, a ready-made matrix, taken basically from SLALIB 
418  _galToJ2000Matrix = threeToSix(numpy.transpose(numpy.array([ 
419          [-0.054875539695716, -0.873437107995315, -0.483834985836994], 
420          [ 0.494109453305607, -0.444829589431879,  0.746982251810510], 
421          [-0.867666135847849, -0.198076386130820,  0.455983795721093]]))) 
422   
423   
424  ############### Supergalactic coordinates 
425   
426  _galToSupergalTrafo = sphermath.computeTransMatrixFromPole( 
427          (47.37*DEG, 6.32*DEG), (137.37*DEG, 0)) 
428  _galToSupergalMatrix = threeToSix(_galToSupergalTrafo) 
429   
430   
431  ############### Ecliptic coordinates 
432   
433 -def _getEclipticMatrix(epoch): 
434          """returns the rotation matrix from equatorial to ecliptic at datetime epoch. 
435   
436          Strictly, epoch should be a TDB. 
437          """ 
438          t = times.getSeconds(epoch-times.dtJ2000)/common.secsPerJCy 
439          obliquity = (84381.448+(-46.8150+(-0.00059+0.001813*t)*t)*t)*ARCSEC 
440          return sphermath.getRotX(obliquity) 
441   
442 -def _getFromEclipticMatrix(fromNode, toNode, sixTrans): 
443          return threeToSix(numpy.transpose(_getEclipticMatrix(fromNode[1]))) 
444   
445 -def _getToEclipticMatrix(fromNode, toNode, sixTrans): 
446          emat = _getEclipticMatrix(fromNode[1]) 
447          return threeToSix(emat) 
448   
449   
450  ############### ICRS a.k.a. Hipparcos 
451  # This is all parallel to IAU sofa, i.e. no zonal corrections, etc. 
452  # From FK5hip 
453   
454 -def cross(vec1, vec2): 
455          """returns the cross product of two 3-vectors. 
456   
457          This should really be somewhere else... 
458          """ 
459          return numpy.array([ 
460                  vec1[1]*vec2[2]-vec1[2]*vec2[1], 
461                  vec1[2]*vec2[0]-vec1[0]*vec2[2], 
462                  vec1[0]*vec2[1]-vec1[1]*vec2[0], 
463          ]) 
464   
465  # Compute transformation from orientation of FK5 
466  _fk5ToICRSMatrix = sphermath.getMatrixFromEulerVector( 
467          numpy.array([-19.9e-3, -9.1e-3, 22.9e-3])*ARCSEC) 
468  _icrsToFK5Matrix = numpy.transpose(_fk5ToICRSMatrix) 
469   
470  # Spin of FK5 in FK5 system 
471  _fk5SpinFK5 = numpy.array([-0.30e-3, 0.60e-3, 0.70e-3])*ARCSEC/365.25 
472  # Spin of FK5 in ICRS 
473  _fk5SpinICRS = numpy.dot(_fk5ToICRSMatrix, _fk5SpinFK5) 
474   
475 -def fk5ToICRS(sixTrans, svFk5): 
476          """returns a 6-vector in ICRS for a 6-vector in FK5 J2000. 
477          """ 
478          spatial = numpy.dot(_fk5ToICRSMatrix, svFk5[:3]) 
479          vel = numpy.dot(_fk5ToICRSMatrix, 
480                  svFk5[3:]+cross(svFk5[:3], _fk5SpinFK5)) 
481          return numpy.concatenate((spatial, vel)) 
482   
483   
484 -def icrsToFK5(sixTrans, svICRS): 
485          """returns a 6-vector in FK5 J2000 for an ICRS 6-vector. 
486          """ 
487          spatial = numpy.dot(_icrsToFK5Matrix, svICRS[:3]) 
488          corrForSpin = svICRS[3:]-cross(svICRS[:3], _fk5SpinICRS) 
489          vel = numpy.dot(_icrsToFK5Matrix, corrForSpin) 
490          return numpy.concatenate((spatial, vel)) 
491   
492   
493  ############### Reference positions 
494  # XXX TODO: We don't transform anything here.  Yet.  This will not 
495  # hurt for moderate accuracy requirements in the stellar and 
496  # extragalactic regime but makes this library basically useless for 
497  # solar system work. 
498   
499 -def _transformRefpos(sixTrans, sixVec): 
500          return sixVec 
501   
502   
503  ############### Top-level code 
504   
505   
506 -def _Constant(val): 
507          """returns a transform factory always returning val. 
508          """ 
509          return lambda fromSTC, toSTC, sixTrans: val 
510   
511   
512  # transforms are triples of fromNode, toNode, transform factory.  Due to 
513  # the greedy nature of your "virtual graph" searching, it's generally a 
514  # good idea to put more specific transforms further up. 
515   
516  _findTransformsPath = _makeFindPath([ 
517          (("FK4", times.dtB1950, SAME), ("FK5", times.dtJ2000, SAME), 
518                  _Constant(fk4ToFK5)), 
519          (("FK5", times.dtJ2000, SAME), ("FK4", times.dtB1950, SAME), 
520                  _Constant(fk5ToFK4)), 
521          (("FK5", times.dtJ2000, SAME), ("GALACTIC_II", ANYVAL, SAME), 
522                  _Constant(la.inv(_galToJ2000Matrix))), 
523          (("GALACTIC_II", ANYVAL, SAME), ("FK5", times.dtJ2000, SAME), 
524                  _Constant(_galToJ2000Matrix)), 
525          (("FK4", times.dtB1950, SAME), ("GALACTIC_II", ANYVAL, SAME), 
526                  _Constant(_b1950ToGalMatrix)), 
527          (("GALACTIC_II", ANYVAL, SAME), ("FK4", times.dtB1950, SAME), 
528                  _Constant(la.inv(_b1950ToGalMatrix))), 
529          (("GALACTIC_II", ANYVAL, SAME), ("SUPER_GALACTIC", ANYVAL, SAME), 
530                  _Constant(_galToSupergalMatrix)), 
531          (("SUPER_GALACTIC", ANYVAL, SAME), ("GALACTIC_II", ANYVAL, SAME), 
532                  _Constant(la.inv(_galToSupergalMatrix))), 
533          (("FK5", ANYVAL, SAME), ("FK5", times.dtJ2000, SAME), 
534                  _getIAU1976PrecMatrix), 
535          (("FK4", ANYVAL, SAME), ("FK4", times.dtB1950, SAME), 
536                  _getNewcombPrecMatrix), 
537          (("FK5", ANYVAL, SAME), ("FK5", ANYVAL, SAME), 
538                  _getIAU1976PrecMatrix), 
539          (("FK4", ANYVAL, SAME), ("FK4", ANYVAL, SAME), 
540                  _getNewcombPrecMatrix), 
541          (("ECLIPTIC", SAME, SAME), ("FK5", SAME, SAME), 
542                  _getFromEclipticMatrix), 
543          (("FK5", SAME, SAME), ("ECLIPTIC", SAME, SAME), 
544                  _getToEclipticMatrix), 
545          (("FK5", times.dtJ2000, SAME), ("ICRS", ANYVAL, SAME), 
546                  _Constant(fk5ToICRS)), 
547          (("ICRS", ANYVAL, SAME), ("FK5", times.dtJ2000, SAME), 
548                  _Constant(icrsToFK5)), 
549          ((SAME, SAME, ANYVAL), (SAME, SAME, ANYVAL), 
550                  _Constant(_transformRefpos)), 
551  ]) 
552   
553   
554  _precessionFuncs = set([_getNewcombPrecMatrix, _getIAU1976PrecMatrix]) 
555   
556 -def _contractPrecessions(toContract): 
557          """contracts the precessions in toContract. 
558           
559          No checks done.  This is only intended as a helper for _simplifyPath. 
560          """ 
561          return toContract[0][0], toContract[-1][1], toContract[0][-1] 
562   
563   
564 -def _simplifyPath(path): 
565          """tries to simplify path by contracting mulitple consecutive precessions. 
566   
567          These come in since our path finding algorithm sucks.  This is mainly 
568          done for numerical reasons since the matrices would be contracted for 
569          computation anyway. 
570          """ 
571  # Sorry about this complex mess.  Maybe we want a more general optimization 
572  # framework. 
573          if path is None: 
574                  return path 
575          newPath, toContract = [], [] 
576          curPrecFunc = None 
577          for t in path: 
578                  if curPrecFunc: 
579                          if t[-1] is curPrecFunc: 
580                                  toContract.append(t) 
581                          else: 
582                                  newPath.append(_contractPrecessions(toContract)) 
583                                  if t[-1] in _precessionFuncs: 
584                                          curPrecFunc, toContract = t[-1], [t] 
585                                  else: 
586                                          curPrecFunc, toContract = None, [] 
587                                          newPath.append(t) 
588                  else: 
589                          if t[-1] in _precessionFuncs: 
590                                  curPrecFunc = t[-1] 
591                                  toContract = [t] 
592                          else: 
593                                  newPath.append(t) 
594          if toContract: 
595                  newPath.append(_contractPrecessions(toContract)) 
596          return newPath 
597   
598 -def _contractMatrices(ops): 
599          """combines consecutive numpy.matrix instances in the sequence 
600          ops by dot-multiplying them. 
601          """ 
602          newSeq, curMat = [], None 
603          for op in ops: 
604                  if isinstance(op, numpy.ndarray): 
605                          if curMat is None: 
606                                  curMat = op 
607                          else: 
608                                  curMat = numpy.dot(curMat, op) 
609                  else: 
610                          if curMat is not None: 
611                                  newSeq.append(curMat) 
612                                  curMat = None 
613                          newSeq.append(op) 
614          if curMat is not None: 
615                  newSeq.append(curMat) 
616          return newSeq 
617           
618   
619 -def _pathToFunction(trafoPath, sixTrans): 
620          """returns a function encapsulating all operations contained in 
621          trafoPath. 
622   
623          The function receives and returns a 6-vector.  trafoPath is altered. 
624          """ 
625          trafoPath.reverse() 
626          steps = _contractMatrices([factory(srcTrip, dstTrip, sixTrans) 
627                  for srcTrip, dstTrip, factory in trafoPath]) 
628          expr = [] 
629          for index, step in enumerate(steps): 
630                  if isinstance(step, numpy.ndarray): 
631                          expr.append("numpy.dot(steps[%d], "%index) 
632                  else: 
633                          expr.append("steps[%d](sixTrans, "%index) 
634          vars = {"steps": steps, "numpy": numpy} 
635          exec ("def transform(sv, sixTrans): return %s"% 
636                  "".join(expr)+"sv"+(")"*len(expr))) in vars 
637          return vars["transform"] 
638   
639   
640 -def nullTransform(sv, sixTrans): 
641          return sv 
642   
643   
644  @memoized 
645 -def getTrafoFunction(srcTriple, dstTriple, sixTrans): 
646          """returns a function that transforms 6-vectors from the system 
647          described by srcTriple to the one described by dstTriple. 
648   
649          The triples consist of (system, equinox, refpoint). 
650   
651          If no transformation function can be produced, the function raises 
652          an STCValueError. 
653   
654          sixTrans is a sphermath.SVConverter instance, used here for communication 
655          of input details and user preferences. 
656          """ 
657          # special case the identity since it's indistingishable from a failed 
658          # search otherwise 
659          if srcTriple==dstTriple: 
660                  return nullTransform 
661          trafoPath = _simplifyPath(_findTransformsPath(srcTriple, dstTriple)) 
662          if trafoPath is None: 
663                  raise common.STCValueError("Cannot find a transform from %s to %s"%( 
664                          srcTriple, dstTriple)) 
665          return _pathToFunction(trafoPath, sixTrans) 
666