gavo.stc.sphermath

1 """ 2 Spherical geometry and 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 import new 13 import numpy 14 15 from gavo.stc import common 16 from gavo.stc import units 17 18 from gavo.utils.mathtricks import cartToSpher, spherToCart #noflake 19 20 # filled in for distances not given, in rad (units also insert this for 21 # parallaxes too small) 22 defaultDistance = units.maxDistance*units.onePc/units.oneAU 23 24 # Units spherical coordinates have to be in for transformation to/from 25 # 6-vectors; this is all done in the bowels of SVConverter 26 _svPosUnit = ("rad", "rad", "AU") 27 _svVPosUnit = ("rad", "rad", "AU") 28 _svVTimeUnit = ("d", "d", "d") 29 # Light speed in AU/d 30 _lightAUd = 86400.0/499.004782

31 32 33 -def getRotX(angle):

34 """returns a 3-rotation matrix for rotating angle radians around x. 35 """ 36 c, s = math.cos(angle), math.sin(angle) 37 return numpy.array([[1, 0, 0], [0, c, s], [0, -s, c]])

38

39 40 -def getRotY(angle):

41 """returns a 3-rotation matrix for rotating angle radians around y. 42 """ 43 c, s = math.cos(angle), math.sin(angle) 44 return numpy.array([[c, 0, -s], [0, 1, 0], [s, 0, c]])

45

46 47 -def getRotZ(angle):

48 """returns a 3-rotation matrix for rotating angle radians around u. 49 """ 50 c, s = math.cos(angle), math.sin(angle) 51 return numpy.array([[c, s, 0], [-s, c, 0], [0, 0, 1]])

52

53 54 -def getMatrixFromEulerAngles(z1, x, z2):

55 """returns a 3-rotation matrix for the z-x-z Euler angles. 56 57 There are some functions to obtain such angles below. 58 """ 59 return numpy.dot( 60 numpy.dot(getRotZ(z2), getRotX(x)), getRotZ(z1))

61

62 63 -def getEulerAnglesFromMatrix(matrix):

64 """returns zxz Euler angles from a rotation matrix. 65 66 This is improvised, and someone should look up a numerically sound way 67 to do this. 68 """ 69 # while I cannot see why, this clearly is broken in some way. 70 z1 = math.atan2(matrix[2,0], -matrix[2,1]) 71 z2 = math.atan2(matrix[0,2], matrix[1,2]) 72 x = math.atan2(math.sqrt(matrix[2,0]**2+matrix[2,1]**2), matrix[2,2]) 73 return z1, x, z2

74

75 76 -def getMatrixFromEulerVector(eulerVector):

77 """returns a rotation matrix for an Euler vector. 78 79 An euler vector gives the rotation axis, its magnitude the angle in rad. 80 81 This function is a rip-off of SOFA's rv2m. 82 83 eulerVector is assumed to be a numpy array. 84 """ 85 x, y, z = eulerVector 86 phi = math.sqrt(x**2+y**2+z**2) 87 sp, cp = math.sin(phi), math.cos(phi) 88 f = 1-cp 89 if phi!=0: 90 x, y, z = eulerVector/phi 91 return numpy.array([ 92 [x**2*f+cp, x*y*f+z*sp, x*z*f-y*sp], 93 [y*x*f-z*sp, y**2*f+cp, y*z*f+x*sp], 94 [z*x*f+y*sp, z*y*f-x*sp, z**2*f+cp]])

95

96 97 -def computeTransMatrixFromPole(poleCoo, longZeroCoo, changeHands=False):

98 """returns a transformation matrix to transform from the reference 99 system into a rotated system. 100 101 The rotated system is defined by its pole, the spherical coordinates 102 at which it has longitude zero and whether or not it is right handed. 103 104 All angles are in rad. 105 """ 106 # when moving from numpy here, replace this with the like-named 107 # function from utils.mathtricks. 108 x = spherToCart(*longZeroCoo) 109 z = spherToCart(*poleCoo) 110 if abs(numpy.dot(x, z))>1e-5: 111 raise common.STCValueError("%s and %s are not valid pole/zero points for" 112 " a rotated coordinate system"%(poleCoo, longZeroCoo)) 113 y = (z[1]*x[2]-z[2]*x[1], z[2]*x[0]-z[0]*x[2], z[0]*x[1]-z[1]*x[0]) 114 if changeHands: 115 y = (-y[0], -y[1], -y[2]) 116 return numpy.array([x,y,z]) 117

118 119 -def vabs(naVec):

120 return math.sqrt(numpy.dot(naVec, naVec))

121

122 123 -def _ensureSphericalFrame(coo):

124 """raises an error if coo's frame is not suitable for holding spherical 125 coordinates. 126 XXX TODO: assert spatial and vel coos have the same frame, etc. 127 """ 128 if not coo.frame.isSpherical(): 129 raise common.STCValueError("%s is not a valid frame for transformable" 130 " spherical coordinates."%(coo.frame))

131

132 133 ############## Code for handling relativistic correction in 6-Vector conversion 134 # Warning: This has massive issues with numerics. Don't use. 135 136 -def _pleaseEinsteinToSpher(sv):

137 """undoes relativistic corrections from 6-vector sv. 138 139 This follows sofa's pvstar. sv is changed in place. 140 """ 141 radialProj, radialV, tangentialV = _decomposeRadial(sv[:3], sv[3:]) 142 betaRadial = radialProj/_lightAUd 143 betaTangential = vabs(tangentialV)/_lightAUd 144 145 d = 1.0+betaRadial 146 w = 1.0-betaRadial**2-betaTangential**2 147 if d==0.0 or w<0: 148 return 149 delta = math.sqrt(w)-1.0 150 if betaRadial==0: 151 radialV = (betaRadial-delta)/(betaRadial*d)*radialV 152 sv[3:] = 1/d*radialV

153

154 155 -def _decomposeRadial(r, rd):

156 """returns the components of rd radial and tangential to r. 157 """ 158 rUnit = r/vabs(r) 159 radialProj = numpy.dot(rUnit, rd) 160 radialVector = radialProj*rUnit 161 tangentialVector = rd-radialVector 162 return radialProj, radialVector, tangentialVector

163

164 165 -def _solveStumpffEquation(betaR, betaT, maxIter=100):

166 """returns the solution of XXX. 167 168 If the solution fails to converge within maxIter iterations, it 169 raises an STCError. 170 """ 171 curEstR, curEstT = betaR, betaT 172 od, odel, odd, oddel = 0, 0, 0, 0 173 for i in range(maxIter): 174 d = 1.+curEstT 175 delta = math.sqrt(1.-curEstR**2-curEstT**2)-1.0 176 curEstR = d*betaR+delta 177 curEstT = d*betaT 178 if i: # check solution so far after at least one iteration 179 dd = abs(d-od) 180 ddel = abs(delta-odel) 181 if dd==odd and ddel==oddel: 182 break 183 odd = dd 184 oddel = ddel 185 od = d 186 odel = delta 187 else: 188 raise common.STCError( 189 "6-vector relativistic correction failed to converge") 190 return curEstR, curEstT, d, delta

191

192 193 -def _pleaseEinsteinFromSpher(sv):

194 """applies relativistic corrections to the 6-vector sv. 195 196 This follows sofa's starpv. sv is changed in place. 197 """ 198 radialProj, radialV, tangentialV = _decomposeRadial(sv[:3], sv[3:]) 199 betaRadial = radialProj/_lightAUd 200 betaTangential = vabs(tangentialV)/_lightAUd 201 202 betaSR, betaST, d, delta = _solveStumpffEquation(betaRadial, betaTangential) 203 # replace old velocity with velocity in inertial system 204 if betaSR!=0: 205 radialV = (d+delta/betaSR)*radialV 206 sv[3:] = radialV+(d*tangentialV)

207

208 ################# End relativistic code. 209 210 -class SVConverter(object):

211 """A container for the conversion from spherical coordinates 212 to 6-Vectors. 213 214 You create one with an example of your data; these are values 215 and units of STC objects, and everything may be None if it's 216 not given. 217 218 The resulting object has methods to6 taking values 219 like the one provided by you and returning a 6-vector, and from6 220 returning a pair of such values. 221 222 Further, the converter has the attributes distGiven, 223 posdGiven, and distdGiven signifying whether these items are 224 expected or valid on return. If posVals is None, no transforms 225 can be computed. 226 227 The relativistic=True constructior argument requests that the 228 transformation be Lorentz-invariant. Do not use that, though, 229 since there are unsolved numerical issues. 230 231 The slaComp=True constructor argument requests that some 232 operations exterior to the construction are done as slalib does 233 them, rather than alternative approaches. 234 """ 235 posGiven = distGiven = posdGiven = distdGiven = True 236 defaultDistance = units.maxDistance*units.onePc/units.oneAU 237

238 - def __init__(self, posVals, posUnit, velVals=None, velSUnit=None, 239 velTUnit=None, relativistic=False, slaComp=False):

240 self.relativistic, self.slaComp = relativistic, slaComp 241 self._determineFeatures(posVals, velVals) 242 self._computeUnitConverters(posUnit, velSUnit, velTUnit) 243 self._makeTo6() 244 self._makeFrom6()

245

246 - def _determineFeatures(self, posVals, velVals):

247 if posVals is None: 248 # No position has been given; that is possible if only "wiggles" 249 # were specified. We'd have to transform those, too, but we 250 # don't for now. 251 self.distGiven = self.posdGiven = False 252 self._velDefault = None 253 return 254 255 if len(posVals)==2: 256 self.distGiven = False 257 if velVals is None: 258 self.posdGiven = False 259 self._velDefault = None 260 else: 261 if len(velVals)==2: 262 self.distdGiven = False 263 self._velDefault = (0,0) 264 else: 265 self._velDefault = (0,0,0)

266

267 - def _computeUnitConverters(self, posUnit, velSUnit, velTUnit):

268 dims = len(posUnit) 269 self.toSVUnitsPos = units.getVectorConverter(posUnit, 270 _svPosUnit[:dims]) 271 self.fromSVUnitsPos = units.getVectorConverter(posUnit, 272 _svPosUnit[:dims], True) 273 if self.posdGiven: 274 dims = len(velSUnit) 275 self.toSVUnitsVel = units.getVelocityConverter(velSUnit, velTUnit, 276 _svVPosUnit[:dims], _svVTimeUnit[:dims]) 277 self.fromSVUnitsVel = units.getVelocityConverter(velSUnit, velTUnit, 278 _svVPosUnit[:dims], _svVTimeUnit[:dims], True)

279

280 - def _makeTo6(self):

281 velDefault = "" 282 if not self.posdGiven: 283 velDefault = "=None" 284 code = ["def to6(self, pos, vel%s):"%velDefault] 285 code.append(" pos = self.toSVUnitsPos(pos)") 286 if self.posdGiven: 287 code.append(" vel = self.toSVUnitsVel(vel)") 288 if not self.distGiven: 289 code.append(" pos = pos+(defaultDistance,)") 290 if self.posdGiven: 291 if not self.distdGiven: 292 code.append(" vel = vel+(0,)") 293 else: 294 code.append(" vel = (0,0,0)") 295 code.append(" (alpha, delta, r), (alphad, deltad, rd) = pos, vel") 296 code.append(" sa, ca = math.sin(alpha), math.cos(alpha)") 297 code.append(" sd, cd = math.sin(delta), math.cos(delta)") 298 code.append(" x, y = r*cd*ca, r*cd*sa") 299 code.append(" w = r*deltad*sd-cd*rd") 300 code.append(" res = numpy.array([x, y, r*sd," 301 " -y*alphad-w*ca, x*alphad-w*sa, r*deltad*cd+sd*rd])") 302 if self.relativistic: 303 code.append(" _pleaseEinsteinFromSpher(res)") 304 code.append(" return res") 305 l = locals() 306 exec "\n".join(code) in globals() , l 307 self.to6 = new.instancemethod(l["to6"], self)

308

309 - def _makeFrom6(self):

310 code = ["def from6(self, sv):"] 311 if self.relativistic: 312 code.append(" _pleaseEinsteinFromSpher(sv)") 313 code.append(" pos, vel = self._svToSpherRaw(sv)") 314 if not self.distGiven: 315 code.append(" pos = pos[:2]") 316 code.append(" pos = self.fromSVUnitsPos(pos)") 317 if self.posdGiven: 318 if not self.distdGiven: 319 code.append(" vel = vel[:2]") 320 code.append(" vel = self.fromSVUnitsVel(vel)") 321 else: 322 code.append(" vel = None") 323 code.append(" return pos, vel") 324 l = locals() 325 exec "\n".join(code) in globals() , l 326 self.from6 = new.instancemethod(l["from6"], self)

327

328 - def _svToSpherRaw(self, sv):

329 """returns spherical position and velocity vectors for the cartesian 330 6-vector sv. 331 332 This is based on SOFA's pv2s. 333 """ 334 x, y, z, xd, yd, zd = sv 335 rInXY2 = x**2+y**2 336 r2 = rInXY2+z**2 337 rw = rTrue = math.sqrt(r2) 338 339 if rTrue==0.: # pos is null: use velocity for position 340 x, y, z = sv[3:] 341 rInXY2 = x**2+y**2 342 r2 = rInXY2+z**2 343 rw = math.sqrt(r2) 344 345 rInXY = math.sqrt(rInXY2) 346 xyp = x*xd+y*yd 347 radialVel = 0 348 if rw!=0: 349 radialVel = xyp/rw+z*zd/rw 350 351 if rInXY2!=0.: 352 theta = math.atan2(y, x) 353 if abs(theta)<1e-12: # null out to avoid wrapping to 2 pi 354 theta = 0 355 if theta<0: 356 theta += 2*math.pi 357 posValues = (theta, math.atan2(z, rInXY), rTrue) 358 velValues = ((x*yd-y*xd)/rInXY2, (zd*rInXY2-z*xyp)/(r2*rInXY), radialVel) 359 else: 360 phi = 0 361 if z!=0: 362 phi = math.atan2(z, rInXY) 363 posValues = (0, phi, rTrue) 364 velValues = (0, 0, radialVel) 365 return posValues, velValues

366

367 - def getPlaceTransformer(self, sixTrafo):

368 """returns a function that transforms 2- or 3-spherical coordinates 369 using the 6-vector transformation sixTrafo. 370 371 Regardless of whether we expect distances, the returned function always 372 returns vectors of the same dimensionality as were passed in. 373 374 This is used when transforming areas. 375 """ 376 def sTrafo(pos): 377 make2 = len(pos)==2 378 if self.distGiven and make2: 379 pos = pos+(defaultDistance,) 380 res, _ = self.from6(sixTrafo(self.to6(pos, self._velDefault), self)) 381 if make2: 382 res = res[:2] 383 return res

384 return sTrafo

385

386 - def getVelocityTransformer(self, sixTrafo, basePos):

387 """returns a function that transforms velocities using sixTrafos 388 at basePos. 389 390 This is used when tranforming velocity intervals. 391 """ 392 def sTrafo(vel): 393 return self.from6(sixTrafo(self.to6(basePos, vel), self))[1]

394 return sTrafo 395 396 @classmethod

397 - def fromSTC(cls, stc, **kwargs):

398 """returns a new 6-vector transform for coordinates in STC. 399 """ 400 pos, vel, posUnit, velUnitS, velUnitT = (None,)*5 401 if stc.place: 402 pos = stc.place.value 403 posUnit = stc.place.unit 404 if pos is None: # if areas are there, fake a position of their 405 if stc.areas: # dimensionality 406 pos = (0,)*stc.place.frame.nDim 407 if stc.velocity: 408 vel = stc.velocity.value 409 velUnitS = stc.velocity.unit 410 velUnitT = stc.velocity.velTimeUnit 411 return cls(pos, posUnit, vel, velUnitS, velUnitT, **kwargs)

412

413 414 -def toSpherical(threeVec):

415 """returns spherical coordinates for a cartesian 3-vector. 416 417 threeVec needs not be normalized. 418 """ 419 rho = math.sqrt(threeVec[0]**2+threeVec[1]**2) 420 long = math.atan2(threeVec[1], threeVec[0]) 421 lat = math.atan2(threeVec[2], rho) 422 return long, lat

423

424 425 -def toThreeVec(long, lat):

426 """returns a cartesian 3-vector for longitude and latitude. 427 """ 428 return (math.cos(lat)*math.cos(long), 429 math.cos(lat)*math.sin(long), 430 math.sin(lat))

431

Source Code for Module gavo.stc.sphermath