The BCC Coordinate System

miller_101.py

"""

a demo of how to convert from cubic to BCC lattice, analog to the tetrahedral
or FCC lattice, but without the headache of having to manage an interlaced grid
structure as with the usual approach, at the expense of a diagonally skewed 
non-cubic bounding box in cartesian coordinates.

an integer conversion from bcc to cartesian coordinates is as easy as

x = -tx+ty+tz
y =  tx-ty+tz
z =  tx+ty-tz

equivalently, the conversion from cartesian to bcc coordinates goes

tx = (y+z)/2
ty = (z+x)/2
tz = (x+y)/2

you may notice that this is exactly the inverse of the fcc transform.

the vertex neighborhood of any coordinate consists of 14 points, and
can be reached by the bcc coordinates

( -1 -1 -1 )
( -1  0  0 )
(  0 -1  0 )
(  0  0 -1 )
(  1  1  1 )
(  1  0  0 )
(  0  1  0 )
(  0  0  1 )

(  0  1  1 )
(  0 -1 -1 )
(  1  0  1 )
( -1  0 -1 )
(  1  1  0 )
( -1 -1  0 )

which in integer cartesian coordinates are the diagonals and the orthogonals
of the BCC lattice:

( -1 -1 -1 )
(  1 -1 -1 )
( -1  1 -1 )
( -1 -1  1 )
(  1  1  1 )
( -1  1  1 )
(  1 -1  1 )
(  1  1 -1 )

(  2  0  0 )
( -2  0  0 )
(  0  2  0 )
(  0 -2  0 )
(  0  0  2 )
(  0  0 -2 )

about distance metrics:

the cartesian distance vector between two points A and B in the lattice is

x = -Bx+By+Bz+Ax-Ay-Az
y =  Bx-By+Bz-Ax+Ay-Az
z =  Bx+By-Bz-Ax-Ay+Az

which can also purely be done in bcc metrics

tx = Bx-Ax
ty = By-Ay
tz = Bz-Az

likewise, scaling uniformity is preserved; the lattice is invariant to
translations and scale, but rotations are not preserved; An axis/angle rotation
by Rodrigues' formula nevertheless works if the bcc dot and cross
products are used (see code).

the euclidean norm (length) of a bcc vector is

|v| = sqrt( (y-z)^2 + (x-z)^2 + (x-y)^2 + x^2 + y^2 + z^2 )

or 

|v| = sqrt( 3*(x^2 + y^2 + z^2) - 2*(x*(y + z) + y*z) )

or

|v| = sqrt( 3*(x + y + z)^2 - 8*(x*(y + z) + y*z) )

which saves 1 addition, 2 subtractions, 1 multiplication in comparison.

the cross product of two bcc vectors in cartesian coordinates is

p = u x v

px = 2*(Ux*(Vy-Vz) + Vx*(Uz-Uy))
py = 2*(Uy*(Vz-Vx) + Vy*(Ux-Uz))
pz = 2*(Uz*(Vx-Vy) + Vz*(Uy-Ux))

remaining in bcc coordinates, the cross product is

px = Ux*(  Vy -   Vz) + Uy*(2*Vz -   Vx) + Uz*(  Vx - 2*Vy)
py = Ux*(  Vy - 2*Vz) + Uy*(  Vz -   Vx) + Uz*(2*Vx -   Vy)
pz = Ux*(2*Vy -   Vz) + Uy*(  Vz - 2*Vx) + Uz*(  Vx -   Vy)

the dot product of two bcc vectors is

q = u . v

q = Ux*(3*Vx-Vy-Vz) + Uy*(3*Vy-Vx-Vz) + Uz*(3*Vz-Vx-Vy)

"""

from math import *

def norm(x,y,z):
    return (x*x+y*y+z*z)**0.5
    
def bccnorm(x,y,z):
    xyz = x+y+z
    return (3*xyz*xyz-8*(x*(z+y)+y*z))**0.5

def bcc2cart(tx,ty,tz):
    return (-tx+ty+tz),(tx-ty+tz),(tx+ty-tz)

def cart2bcc(x,y,z):
    return (y+z)*0.5,(x+z)*0.5,(x+y)*0.5

coords = [
( -1,-1,-1 ),
( -1, 0, 0 ),
(  0,-1, 0 ),
(  0, 0,-1 ),
(  1, 1, 1 ),
(  1, 0, 0 ),
(  0, 1, 0 ),
(  0, 0, 1 ),

(  0, 1, 1 ),
(  0,-1,-1 ),
(  1, 0, 1 ),
( -1, 0,-1 ),
(  1, 1, 0 ),
( -1,-1, 0 ),
]

for x,y,z in coords:
    print "( {:2d} {:2d} {:2d} )".format(x,y,z)
print
for x,y,z in coords:
    print bccnorm(x,y,z)
    print bcc2cart(x,y,z)

U = (1,2,3)
V = (0,4,1)

def cross(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    return Uy*Vz - Uz*Vy, Uz*Vx - Ux*Vz, Ux*Vy - Uy*Vx

def dot(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    return Ux*Vx+Uy*Vy+Uz*Vz

def bccdot(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    return Ux*(3*Vx-Vy-Vz) + Uy*(3*Vy-Vx-Vz) + Uz*(3*Vz-Vx-Vy)

def bcccross2cart(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V
    px = 2*(Ux*(Vy-Vz) + Vx*(Uz-Uy))
    py = 2*(Uy*(Vz-Vx) + Vy*(Ux-Uz))
    pz = 2*(Uz*(Vx-Vy) + Vz*(Uy-Ux))
    return px,py,pz

def bcccross(U,V):
    Ux,Uy,Uz = U
    Vx,Vy,Vz = V    
    px = Ux*(  Vy -   Vz) + Uy*(2*Vz -   Vx) + Uz*(  Vx - 2*Vy)
    py = Ux*(  Vy - 2*Vz) + Uy*(  Vz -   Vx) + Uz*(2*Vx -   Vy)
    pz = Ux*(2*Vy -   Vz) + Uy*(  Vz - 2*Vx) + Uz*(  Vx -   Vy)
    return px,py,pz
    
def taxicab(x,y,z):
    return abs(x)+abs(y)+abs(z)

def bccaxisangle(a,U,W):
    s,c = sin(a),cos(a)
    Sx,Sy,Sz = bcccross(W,U)
    t = (1-c)*bccdot(W,U)
    Ux,Uy,Uz = U
    Wx,Wy,Wz = W
    vx = c*Ux + s*Sx + t*Wx
    vy = c*Uy + s*Sy + t*Wy
    vz = c*Uz + s*Sz + t*Wz
    return vx,vy,vz

def axisangle(a,U,W):
    s,c = sin(a),cos(a)
    Sx,Sy,Sz = cross(W,U)
    t = (1-c)*dot(W,U)
    Ux,Uy,Uz = U
    Wx,Wy,Wz = W
    vx = c*Ux + s*Sx + t*Wx
    vy = c*Uy + s*Sy + t*Wy
    vz = c*Uz + s*Sz + t*Wz
    return vx,vy,vz

print(dot((1,2,3),(6,5,4)))
print(dot((1,2,3),(4,5,6)))
print(bccdot(cart2bcc(1,2,3),cart2bcc(6,5,4)))
print(bccdot(cart2bcc(1,2,3),cart2bcc(4,5,6)))
print(cross((1,2,3),(6,5,4)))
print(bcccross2cart(cart2bcc(1,2,3),cart2bcc(6,5,4)))
print(bcc2cart(*bcccross(cart2bcc(1,2,3),cart2bcc(6,5,4))))

print(norm(4,5,6))
print(norm(*axisangle(radians(12.3),(4,5,6),(0,-0.70710678,0.70710678))))
print(bccnorm(*bccaxisangle(radians(12.3),cart2bcc(4,5,6),cart2bcc(0,-0.70710678,0.70710678))))