import sys
import subprocess
import math as mth
import cmath as cmt
import numpy as npy
import scipy.linalg as sla

dt = 0.08
iterations = 100
state_num = 1

H = npy.array([[1,0],[0,2]])

if len(sys.argv) == 3:
    alpha = complex(sys.argv[1])
    beta = complex(sys.argv[2])

def init_state(i):
    alpha = 1
    beta  = 1
    if state_num > 1:
        alpha = i / state_num
        beta  = (state_num - i) / state_num
    return (alpha, beta)

def time_evolution(state, dt = dt):
    return npy.dot(state, sla.expm(-1j * dt * H))

def bloch_map(state):
    temp  = npy.conj(state[0] / state[1])
    theta = 2 * npy.arctan(abs(temp))
    phi = cmt.phase(temp)

    return (phi, theta)

def sphere2cart(phi, theta):
        return [
             mth.sin(theta) * mth.cos(phi),
             mth.sin(theta) * mth.sin(phi),
             -1 * mth.cos(theta)
        ]

f = open("data", "w")

states = []
for i in range(state_num):
    (alpha, beta) = init_state(i)
    norm  = npy.linalg.norm([alpha, beta])
    state = npy.array([alpha / norm, beta / norm])
    states.append(state)


for i in range(iterations):
    for j in range(state_num):
        (phi, theta) = bloch_map(states[j])
        coords = sphere2cart(phi, theta)
        colour = i / iterations
        if state_num > 1:
            colour = j / state_num
        f.write(f"{coords[0]}; {coords[1]}; {coords[2]}; {colour}\n")
        states[j] = time_evolution(states[j])

f.close()

import os

path = 'gnuplot/'
if os.name == 'nt':
    path.replace('/','\\')
subprocess.run(["gnuplot",
    "-e", f"states={state_num};",
    "-e", f"start={iterations};" ,
    "-e", f"end={iterations};",
    "-c", path + "animation.plt"])