# Timothy Keller
# S = 1/2, I = 1/2
# Spin 1/2 electron coupled to spin 1/2 nuclei

import numpy as np
from scipy.linalg import expm
from matplotlib.pylab import *
from matplotlib import cm

sigma_x = 0.5*np.r_[[[0, 1],[1, 0]]]
sigma_y = 0.5*np.r_[[[0,-1j],[1j, 0]]]
sigma_z = 0.5*np.r_[[[1, 0],[0, -1]]]
Identity = np.eye(2)


Sx = np.kron(sigma_x, Identity)
Sy = np.kron(sigma_y, Identity)
Sz = np.kron(sigma_z, Identity)

Ix = np.kron(Identity, sigma_x)
Iy = np.kron(Identity, sigma_y)
Iz = np.kron(Identity, sigma_z)

SxIx = np.kron(sigma_x,sigma_z)

SxIx2 = np.dot(Sx,Iz)

print(SxIx)
print(SxIx2)
print(np.allclose(SxIx,SxIx2))
omega_S = 1.76e11 # rad / (s * T)
omega_I = 267.522e6 # rad / (s * T)
Aiso = 2*np.pi * 50.e6 # Isotropic Hyperfine coupling rad / s

B0 = 0.35# T

H = omega_S/(2.*np.pi)*B0*Sz + omega_I/(2.*np.pi)*B0*Iz + Aiso * np.dot(Sz,Iz)
#H = omega_S/(2.*np.pi)*B0*Sz + omega_I/(2.*np.pi)*B0*Iz + Aiso * (np.dot(Sx,Ix) + np.dot(Sy,Iy) + np.dot(Sz,Iz))

print('Hamiltonian:')
print(H)
out = np.linalg.eig(H)

E = out[0]
print(E)

E12 = E[0] - E[1]
E34 = E[2] - E[3]
E13 = E[0] - E[2]
E24 = E[1] - E[3]
print(E12)
print(E34)
print(E13)
print(E24)

print('Nuclear')
print('%0.05f MHz'%(E12 / 1e6))
print('%0.05f MHz'%(E34 / 1e6))
print('Electron')
print('%0.05f GHz'%(E13 / 1e9))
print('%0.05f GHz'%(E24 / 1e9))

matshow(abs(H), cmap = cm.jet)
title('Hamiltonian')
show()

[((13, 11, 13, 20), 'numpy.eye', 'np.eye', ({(13, 18, 13, 19): '2'}, {}), '(2)', True, 'import numpy as np\n'), ((16, 5, 16, 31), 'numpy.kron', 'np.kron', ({(16, 13, 16, 20): 'sigma_x', (16, 22, 16, 30): 'Identity'}, {}), '(sigma_x, Identity)', True, 'import numpy as np\n'), ((17, 5, 17, 31), 'numpy.kron', 'np.kron', ({(17, 13, 17, 20): 'sigma_y', (17, 22, 17, 30): 'Identity'}, {}), '(sigma_y, Identity)', True, 'import numpy as np\n'), ((18, 5, 18, 31), 'numpy.kron', 'np.kron', ({(18, 13, 18, 20): 'sigma_z', (18, 22, 18, 30): 'Identity'}, {}), '(sigma_z, Identity)', True, 'import numpy as np\n'), ((20, 5, 20, 31), 'numpy.kron', 'np.kron', ({(20, 13, 20, 21): 'Identity', (20, 23, 20, 30): 'sigma_x'}, {}), '(Identity, sigma_x)', True, 'import numpy as np\n'), ((21, 5, 21, 31), 'numpy.kron', 'np.kron', ({(21, 13, 21, 21): 'Identity', (21, 23, 21, 30): 'sigma_y'}, {}), '(Identity, sigma_y)', True, 'import numpy as np\n'), ((22, 5, 22, 31), 'numpy.kron', 'np.kron', ({(22, 13, 22, 21): 'Identity', (22, 23, 22, 30): 'sigma_z'}, {}), '(Identity, sigma_z)', True, 'import numpy as np\n'), ((24, 7, 24, 31), 'numpy.kron', 'np.kron', ({(24, 15, 24, 22): 'sigma_x', (24, 23, 24, 30): 'sigma_z'}, {}), '(sigma_x, sigma_z)', True, 'import numpy as np\n'), ((26, 8, 26, 21), 'numpy.dot', 'np.dot', ({(26, 15, 26, 17): 'Sx', (26, 18, 26, 20): 'Iz'}, {}), '(Sx, Iz)', True, 'import numpy as np\n'), ((42, 6, 42, 22), 'numpy.linalg.eig', 'np.linalg.eig', ({(42, 20, 42, 21): 'H'}, {}), '(H)', True, 'import numpy as np\n'), ((30, 6, 30, 29), 'numpy.allclose', 'np.allclose', ({(30, 18, 30, 22): 'SxIx', (30, 23, 30, 28): 'SxIx2'}, {}), '(SxIx, SxIx2)', True, 'import numpy as np\n'), ((37, 65, 37, 78), 'numpy.dot', 'np.dot', ({(37, 72, 37, 74): 'Sz', (37, 75, 37, 77): 'Iz'}, {}), '(Sz, Iz)', True, 'import numpy as np\n')]