In [3]:
"""
Author: Wells Bishop
Date: March 14th, 2016
Lorentz Transformation Matrix Calculator

Linear transformation from one reference frame in R^4 to another in R^4.
X = (t, x, y, z) is the 4-vector giving us the time and position of the moving object.
X'= (t', x', y', z') gives us the vector of how the object appears to an observer.
V is the velocity vector of the object with respect to the observer and will greatly
influence the difference between X and X'. As V gets larger and approaches the speed of light,
the matrix A stretches and contracts various aspects of the source to maintain a constant C.

Lengths get shorter in the direction of travel, at least to an outside observr. Things moving at normal speeds 
don't experience large values, but once they start having velocities towards C=1, objects shrink a great deal
and experience time differently then us.

I followed the matrix on this website (http://www.physicspages.com/2011/06/22/lorentz-transformations-in-three-dimensions/)
and used the knowledge I gained from my physics 351, special relativity to write this program.

Input: 
        Position 4-vector X = (t, x, y, z) in objects reference frame
        Velocity 3-vector V = (vx, vy, vz) between object and observer
Output:
        X' = (t', x', y', z') in the observer's reference frame. This is how the obsever, observes the objects 
                            time and position, since the object experiences time differently when the two 
                            have very large velocities differences

"""

import numpy as np
from scipy import linalg

Position = np.zeros(4) #t, x, y, z
Velocity = np.array([1.0, 0.0, 0.0])
C = 3.00*10**8 #speed of light


def Gamma(num):
    """
    Calculates Gamma factor from the velocity of the object with
    respect to a stationary observer.
    Input: Velocity or velocity components (int or array)
    Returns: Gamma
    """
    if type(num) is int:
        v = num
    else:
        v = 0.0
        for el in num:
            v += el**2
        v = np.sqrt(v)

    G = np.sqrt(1-(v/C)**2)
    return 1.0 / G

def Lorentz(Gamma, V):
    """
    Creates a 4x4 Lorentz transformation matrix. Which is a linear transformation
    from R^4 to R^4. (t, x, y, z) to (t', x', y', z')
    """
    A = np.zeros((4,4))#4x4 matrix
    vx = V[0]
    vy = V[1]
    vz = V[2]
    v = np.sqrt(vx**2 + vy**2 + vz**2)
    
    A[0,0] = Gamma
    A[1,1] = 1 + (Gamma-1)*((vx/v)**2)
    A[2,2] = 1 + (Gamma-1)*((vy/v)**2)
    A[3,3] = 1 + (Gamma-1)*((vz/v)**2)
    A[0,1] = A[1,0] = -vx*Gamma
    A[0,2] = A[2,0] = -vy*Gamma
    A[0,3] = A[3,0] = -vz*Gamma
    A[1,2] = A[2,1] = (Gamma-1)*(vx*vy)/(v**2)
    A[1,3] = A[3,1] = (Gamma-1)*(vx*vz)/(v**2)
    A[2,3] = A[3,2] = (Gamma-1)*(vy*vz)/(v**2)
    
    E_val, E_vec = linalg.eigh(A)
    return A , E_val, E_vec
    
    
##Evaluate for a few examples    

def Evaluate(A, X, V):
    """
    Given a matrix A, vector X, solve AX=x
    """
    B = np.zeros(4)
    for j in range(0,4):
        k = 0.0
        for i in range(0,4):  # U = {u1,u2,u3,u4}
            u = A[j,i]*X[i]   # u1 = SUM(A[i]*X[i]) for ith element along the jth row
            k = k + u         # linear combination right here
        B[j] = k
    return B  

X = [0, 1, 0, 0]
V = [0.01, 0, 0]
A, Eval, Evec = Lorentz(Gamma(V), V) 
x = Evaluate(A, X, V) ## Ax = X

print "Examples:"
print "First: In one reference frame the object is 1m long, but now in a new frame of reference at 1% the speed of light: "
print "Original time: {}, New time: {}".format(X[0], x[0])
print "Original x-position: {}, new: {}".format(X[1], x[1])
print "Eigenvalues of Lorentz matrix: {}".format(str(Eval))
print "\n"

X = [0, 1, 0, 0]
V = [0.1, 0, 0]
A, Eval, Evec = Lorentz(Gamma(V), V)
x = Evaluate(A, X, V) ## Ax = X
print "Second: In one reference frame the object flys 1 seconds and is 1m long, but now in a new frame of reference at 10% the speed of light: "
print "Original time: {}, New time: {}".format(X[0], x[0])
print "Original x-position: {}, new: {}".format(X[1], x[1])
print "\n"
print "Since the Lorentz Transformation matrix is a symmetric matrix, we can orthogonally diagonalize it."
##Orthogonaly Diagonalize A
print "The matrix A is a square symmetric matrix and therefore, can be orthogoanlly diagonalized by A = P*D*P(transpose) "
print "where D is a diagonal matrix with eigenvalues along the diagonal, and P has the corresponding eigenvectors as its columns:"

P = Evec
D = np.array([ [Eval[0],0,0,0], [0,Eval[1],0,0], [0,0,Eval[2],0], [0,0,0,Eval[3]] ])
PT = np.transpose(P)


print "A= "
print A
print "D= "
print D
print "P= "
print P


X = [1.0, 1/np.sqrt(2), 1/np.sqrt(2), 0.0] #45 degree angle in the x-y plane
V = [0.0, 0.5, 0.0] 
A, Eval, Evec = Lorentz(Gamma(V), V)
x = Evaluate(A, X, V) ## Ax = X
print "\n"
print "Third: In one reference frame the object flys 1 seconds and is 1m long, 10% speed of light y direction: "
print "Original x-position: {}, new: {}".format(X[1], x[1])
print "Original y-position: {}, new: {}".format(X[2], x[2])
print "Eigenvalues of Lorentz matrix: {}".format(str(Eval))
print "\n"

X = [1.0, 1.0, 0.0, 0.0]
V = [0.995, 0.0, 0.0]
A, Eval, Evec = Lorentz(Gamma(V), V)
x = Evaluate(A, X, V) ## Ax = X
print "Fourth: In one reference frame the object flys 1 seconds and is 1m long, but now in a new frame of reference at 99.5% the speed of light: "
print "Original time: {}, New time: {}".format(X[0], x[0])
print "Original x-position: {}, new: {}".format(X[1], x[1])
print "Eigenvalues of Lorentz matrix: {}".format(str(Eval))
##Orthogonaly Diagonalize A
print "The matrix A is a square symmetric matrix and therefore, can be orthogoanlly diagonalized by A = P*D*P(transpose) "
print "where D is a diagonal matrix with eigenvalues along the diagonal, and P has the corresponding eigenvectors as its columns:"

P = Evec
D = np.array([ [Eval[0],0,0,0], [0,Eval[1],0,0], [0,0,Eval[2],0], [0,0,0,Eval[3]] ])
PT = np.transpose(P)


print "A= "
print A
print "D= "
print D
print "P= "
print P
print "\n"

X = [1.0, 1/np.sqrt(3), 1/np.sqrt(3), 1/np.sqrt(3)]
V = [0.94, 0.0, 0.5]
A, Eval, Evec = Lorentz(Gamma(V), V)
x = Evaluate(A, X, V) ## Ax = X
print "Fifth: In one reference frame the object flys 1 seconds and is a 1m side length cube, but now in a ref frame of 0.95*(The speed of light) in the x-z direction: "
print "Original time: {}, New time: {}".format(X[0], x[0])
print "Original x-position: {}, new: {}".format(X[1], x[1])
print "Original y-position: {}, new: {}".format(X[2], x[2])
print "Original z-position: {}, new: {}".format(X[3], x[3])
print "Eigenvalues of Lorentz matrix: {}".format(str(Eval))
print "Now the cube's looks more like a large flat rectangle but is length contracted into a cube from the observers point of view."
print "\n"

##Orthogonaly Diagonalize A
print "\n"
print "The matrix A is a square symmetric matrix and therefore, can be orthogoanlly diagonalized by A = P*D*P(transpose) "
print "where D is a diagonal matrix with eigenvalues along the diagonal, and P has the corresponding eigenvectors as its columns:"

P = Evec
D = np.array([ [Eval[0],0,0,0], [0,Eval[1],0,0], [0,0,Eval[2],0], [0,0,0,Eval[3]] ])
PT = np.transpose(P)
print "D= "
print D
print "P= "
print P


Examples:
First: In one reference frame the object is 1m long, but now in a new frame of reference at 1% the speed of light: 
Original time: 0, New time: -0.01
Original x-position: 1, new: 1.0
Eigenvalues of Lorentz matrix: [ 0.99  1.    1.    1.01]


Second: In one reference frame the object flys 1 seconds and is 1m long, but now in a new frame of reference at 10% the speed of light: 
Original time: 0, New time: -0.1
Original x-position: 1, new: 1.0


Since the Lorentz Transformation matrix is a symmetric matrix, we can orthogonally diagonalize it.
The matrix A is a square symmetric matrix and therefore, can be orthogoanlly diagonalized by A = P*D*P(transpose) 
where D is a diagonal matrix with eigenvalues along the diagonal, and P has the corresponding eigenvectors as its columns:
A= 
[[ 1.  -0.1  0.   0. ]
 [-0.1  1.   0.   0. ]
 [ 0.   0.   1.   0. ]
 [ 0.   0.   0.   1. ]]
D= 
[[ 0.9  0.   0.   0. ]
 [ 0.   1.   0.   0. ]
 [ 0.   0.   1.   0. ]
 [ 0.   0.   0.   1.1]]
P= 
[[ 0.70710678  0.          0.          0.70710678]
 [ 0.70710678  0.          0.         -0.70710678]
 [ 0.          1.          0.          0.        ]
 [ 0.          0.          1.          0.        ]]


Third: In one reference frame the object flys 1 seconds and is 1m long, 10% speed of light y direction: 
Original x-position: 0.707106781187, new: 0.707106781187
Original y-position: 0.707106781187, new: 0.207106781187
Eigenvalues of Lorentz matrix: [ 0.5  1.   1.   1.5]


Fourth: In one reference frame the object flys 1 seconds and is 1m long, but now in a new frame of reference at 99.5% the speed of light: 
Original time: 1.0, New time: 0.005
Original x-position: 1.0, new: 0.005
Eigenvalues of Lorentz matrix: [ 0.005  1.     1.     1.995]
The matrix A is a square symmetric matrix and therefore, can be orthogoanlly diagonalized by A = P*D*P(transpose) 
where D is a diagonal matrix with eigenvalues along the diagonal, and P has the corresponding eigenvectors as its columns:
A= 
[[ 1.    -0.995 -0.    -0.   ]
 [-0.995  1.     0.     0.   ]
 [-0.     0.     1.     0.   ]
 [-0.     0.     0.     1.   ]]
D= 
[[ 0.005  0.     0.     0.   ]
 [ 0.     1.     0.     0.   ]
 [ 0.     0.     1.     0.   ]
 [ 0.     0.     0.     1.995]]
P= 
[[ 0.70710678  0.          0.          0.70710678]
 [ 0.70710678  0.          0.         -0.70710678]
 [ 0.          1.          0.          0.        ]
 [ 0.          0.          1.          0.        ]]


Fifth: In one reference frame the object flys 1 seconds and is a 1m side length cube, but now in a ref frame of 0.95*(The speed of light) in the x-z direction: 
Original time: 1.0, New time: 0.168615612367
Original x-position: 0.57735026919, new: -0.36264973081
Original y-position: 0.57735026919, new: 0.57735026919
Original z-position: 0.57735026919, new: 0.0773502691896
Eigenvalues of Lorentz matrix: [-0.06470653  1.          1.          2.06470653]
Now the cube's looks more like a large flat rectangle but is length contracted into a cube from the observers point of view.




The matrix A is a square symmetric matrix and therefore, can be orthogoanlly diagonalized by A = P*D*P(transpose) 
where D is a diagonal matrix with eigenvalues along the diagonal, and P has the corresponding eigenvectors as its columns:
D= 
[[-0.06470653  0.          0.          0.        ]
 [ 0.          1.          0.          0.        ]
 [ 0.          0.          1.          0.        ]
 [ 0.          0.          0.          2.06470653]]
P= 
[[-0.70710678  0.          0.          0.70710678]
 [-0.62428505  0.         -0.46961297 -0.62428505]
 [ 0.          1.          0.          0.        ]
 [-0.33206652  0.          0.88287239 -0.33206652]]

In [ ]: