agd.Sphere

This module provides basic conversion utilities to manipulate low-dimensional spheres, and related objects : rotations, quaternions, Pauli matrices, etc

View Source

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4"""
  5This module provides basic conversion utilities to manipulate low-dimensional spheres, 
  6and related objects : rotations, quaternions, Pauli matrices, etc
  7"""
  8
  9from . import AutomaticDifferentiation as ad
 10from . import LinearParallel as lp
 11
 12import numpy as np
 13
 14# -----------------------
 15# Equatorial projection
 16# -----------------------
 17
 18def sphere_from_plane(e):
 19    """Produces a point in the unit sphere by projecting a point in the equator plane."""
 20    e = ad.asarray(e)
 21    e2 = lp.dot_VV(e,e)
 22    return ad.array([1.-e2,*(2*e)])/(1.+e2)
 23def plane_from_sphere(q):
 24    """Produces a point in the equator plane from a point in the unit sphere."""
 25    e2 = (1-q[0])/(1+q[0])
 26    return q[1:]*((1+e2)/2.)
 27
 28# ------------------------
 29# Rotations, dimension 3
 30# ------------------------
 31
 32# Two or three dimensional rotation, defined by angle, and axis in dimension three.
 33rotation = lp.rotation
 34
 35def rotation3_from_sphere3(q):
 36    """Produces the rotation associated with a unit quaternion."""
 37    qr,qi,qj,qk = q
 38    return 2*ad.array([
 39        [0.5-(qj**2+qk**2), qi*qj-qk*qr, qi*qk+qj*qr],
 40        [qi*qj+qk*qr, 0.5-(qi**2+qk**2), qj*qk-qi*qr],
 41        [qi*qk-qj*qr, qj*qk+qi*qr, 0.5-(qi**2+qj**2)]])
 42
 43def sphere3_from_rotation3(r):
 44    """Produces the unit quaternion, with positive real part, associated with a rotation."""
 45    qr = np.sqrt(lp.trace(r)+1.)/2.
 46    qi = (r[2,1]-r[1,2])/(4*qr)
 47    qj = (r[0,2]-r[2,0])/(4*qr)
 48    qk = (r[1,0]-r[0,1])/(4*qr)
 49    return ad.array([qr,qi,qj,qk])
 50
 51def ball3_from_rotation3(r,qRef=None):
 52    """Produces an euclidean point from a rotation, 
 53    selecting in the intermediate step the quaternion 
 54    in the same hemisphere as qRef. (Defaults to southern.)"""
 55    q = sphere3_from_rotation3(r)
 56    if qRef is not None: q[:,lp.dot_VV(q,qRef)<0] *= -1
 57    return plane_from_sphere(q)
 58
 59def rotation3_from_ball3(e): 
 60    """Produces a rotation from an euclidean point. 
 61    Also returns the intermediate quaternion."""
 62    q = sphere_from_plane(e)
 63    return rotation3_from_sphere3(q),q
 64
 65# -----------------------
 66# Rotations, dimension 2
 67# -----------------------
 68
 69# For now, this is skipped, since it is not very useful, and not very consistent 
 70# with the three dimensional case (which involves a double cover). 
 71
 72# def rotation_from_sphere(q):
 73# 	"""
 74# 	- (dimension 3) Produces the rotation associated with a unit quaternion.
 75# 	- (dimension 2) Produces the rotation whose first column is the given unit vector.
 76# 	"""
 77# 	if   np.ndim(q)==2: return rotation2_from_sphere1(q)
 78# 	elif np.ndim(q)==4: return rotation3_from_sphere3(q)
 79# 	else: raise ValueError("Unsupported dimension")
 80
 81def rotation2_from_sphere1(q):
 82	"""Produces the rotation whose first column is the given unit vector."""
 83	c,s = q
 84	return ad.array([[c,-s],[s,c]])
 85
 86def sphere1_from_rotation2(r):
 87	"""Produces the unit vector which is the first column of the given rotation"""
 88	return ad.array([r[0,0],r[1,0]])
 89
 90# -----------------------
 91# Pauli matrices
 92# -----------------------
 93
 94def pauli(a,b,c,d=None):
 95	"""
 96	Pauli matrix. Symmetric if d is None, Hermitian otherwise.
 97	Determinant is $a^2-b^2-c^2-d^2$
 98	"""
 99	if d is None: return ad.array([[a+b,c],[c,a-b]])
100	else: return ad.array([[a+b,c+1j*d],[c-1j*d,a-b]])

def sphere_from_plane(e): View Source

19def sphere_from_plane(e):
20    """Produces a point in the unit sphere by projecting a point in the equator plane."""
21    e = ad.asarray(e)
22    e2 = lp.dot_VV(e,e)
23    return ad.array([1.-e2,*(2*e)])/(1.+e2)

Produces a point in the unit sphere by projecting a point in the equator plane.

def plane_from_sphere(q): View Source

24def plane_from_sphere(q):
25    """Produces a point in the equator plane from a point in the unit sphere."""
26    e2 = (1-q[0])/(1+q[0])
27    return q[1:]*((1+e2)/2.)

Produces a point in the equator plane from a point in the unit sphere.

def rotation(theta, axis=None): View Source

25def rotation(theta,axis=None):
26	"""
27	Dimension 2 : by a given angle.
28	Dimension 3 : by a given angle, along a given axis.
29	Three dimensional rotation matrix, with given axis and angle.
30	Adapted from https://stackoverflow.com/a/6802723
31	"""
32	if axis is None:
33		c,s=np.cos(theta),np.sin(theta)
34		return ad.asarray([[c,-s],[s,c]])
35	else:
36		theta,axis = (ad.asarray(e) for e in (theta,axis))
37		axis = axis / np.linalg.norm(axis,axis=0)
38		theta,axis=fd.common_field((theta,axis),(0,1))
39		a = np.cos(theta / 2.0)
40		b, c, d = -axis * np.sin(theta / 2.0)
41		aa, bb, cc, dd = a * a, b * b, c * c, d * d
42		bc, ad_, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
43		return ad.asarray([
44			[aa + bb - cc - dd, 2 * (bc + ad_), 2 * (bd - ac)],
45			[2 * (bc - ad_), aa + cc - bb - dd, 2 * (cd + ab)],
46			[2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

Dimension 2 : by a given angle. Dimension 3 : by a given angle, along a given axis. Three dimensional rotation matrix, with given axis and angle. Adapted from https://stackoverflow.com/a/6802723

def rotation3_from_sphere3(q): View Source

36def rotation3_from_sphere3(q):
37    """Produces the rotation associated with a unit quaternion."""
38    qr,qi,qj,qk = q
39    return 2*ad.array([
40        [0.5-(qj**2+qk**2), qi*qj-qk*qr, qi*qk+qj*qr],
41        [qi*qj+qk*qr, 0.5-(qi**2+qk**2), qj*qk-qi*qr],
42        [qi*qk-qj*qr, qj*qk+qi*qr, 0.5-(qi**2+qj**2)]])

Produces the rotation associated with a unit quaternion.

def sphere3_from_rotation3(r): View Source

44def sphere3_from_rotation3(r):
45    """Produces the unit quaternion, with positive real part, associated with a rotation."""
46    qr = np.sqrt(lp.trace(r)+1.)/2.
47    qi = (r[2,1]-r[1,2])/(4*qr)
48    qj = (r[0,2]-r[2,0])/(4*qr)
49    qk = (r[1,0]-r[0,1])/(4*qr)
50    return ad.array([qr,qi,qj,qk])

Produces the unit quaternion, with positive real part, associated with a rotation.

def ball3_from_rotation3(r, qRef=None): View Source

52def ball3_from_rotation3(r,qRef=None):
53    """Produces an euclidean point from a rotation, 
54    selecting in the intermediate step the quaternion 
55    in the same hemisphere as qRef. (Defaults to southern.)"""
56    q = sphere3_from_rotation3(r)
57    if qRef is not None: q[:,lp.dot_VV(q,qRef)<0] *= -1
58    return plane_from_sphere(q)

Produces an euclidean point from a rotation, selecting in the intermediate step the quaternion in the same hemisphere as qRef. (Defaults to southern.)

def rotation3_from_ball3(e): View Source

60def rotation3_from_ball3(e): 
61    """Produces a rotation from an euclidean point. 
62    Also returns the intermediate quaternion."""
63    q = sphere_from_plane(e)
64    return rotation3_from_sphere3(q),q

Produces a rotation from an euclidean point. Also returns the intermediate quaternion.

def rotation2_from_sphere1(q): View Source

82def rotation2_from_sphere1(q):
83	"""Produces the rotation whose first column is the given unit vector."""
84	c,s = q
85	return ad.array([[c,-s],[s,c]])

Produces the rotation whose first column is the given unit vector.

def sphere1_from_rotation2(r): View Source

87def sphere1_from_rotation2(r):
88	"""Produces the unit vector which is the first column of the given rotation"""
89	return ad.array([r[0,0],r[1,0]])

Produces the unit vector which is the first column of the given rotation

def pauli(a, b, c, d=None): View Source

 95def pauli(a,b,c,d=None):
 96	"""
 97	Pauli matrix. Symmetric if d is None, Hermitian otherwise.
 98	Determinant is $a^2-b^2-c^2-d^2$
 99	"""
100	if d is None: return ad.array([[a+b,c],[c,a-b]])
101	else: return ad.array([[a+b,c+1j*d],[c-1j*d,a-b]])

Pauli matrix. Symmetric if d is None, Hermitian otherwise. Determinant is $a^2-b^2-c^2-d^2$