import numpy as np
from numpy import pi, r_
import matplotlib.pyplot as plt
from scipy import optimize
import csv
from scipy.optimize import curve_fit
from math import cos,sin,sqrt

class MyData:
    spin=[]
    energy=[]
    tag=[]
mydata=MyData()

with open('tsd1.csv','r') as csv_file:
    file=csv.reader(csv_file)
    for i in file:
        x=i[0]
        mydata.spin.append(float(x))
        y=i[1]
        mydata.energy.append(float(y))
        z=1
        mydata.tag.append(int(z))

#pre-requisites
def moi(angle,k):
    pi=3.141592654
    beta=0.38
    result=1.0/(1+pow(5.0/(16.0*pi),0.5)*beta)*(1-pow(5.0/(4*pi),0.5)*beta*cos(angle*pi/180.0+(2.0/3.0)*pi*k))
    return result
def a1(angle):
    return 1.0/(2.0*moi(angle,1))
def a2(angle):
    return 1.0/(2.0*moi(angle,2))
def a3(angle):
    return 1.0/(2.0*moi(angle,3))
def n1(spin,angle):
    rad3=sqrt(3.0)
    gamma=angle*pi/180.0
    result=(2*spin-1)/(spin*(spin+1))*rad3*(rad3*cos(gamma)+sin(gamma))
    return result
def n2(spin,angle):
    rad3=sqrt(3.0)
    gamma=angle*pi/180.0
    result=(2*spin-1)/(spin*(spin+1))*2*rad3*sin(gamma)
    return result
def t1(spin,angle,singleparticle,izero):
    result=(2*spin-1)*(1.0/izero*a3(angle)-1.0/izero*a1(angle))+2*singleparticle*1/izero*a1(angle)
    return result
def t2(spin,angle,singleparticle,izero):
    result=(2*spin-1)*(1.0/izero*a2(angle)-1.0/izero*a1(angle))+2*singleparticle*1/izero*a1(angle)
    return result
def t3(spin,angle,singleparticle,izero,v):
    result=(2*singleparticle-1)*(1.0/izero*a3(angle)-1.0/izero*a1(angle))+2.0*spin*1.0/izero*a1(angle)+v*n1(singleparticle,angle)
    return result
def t4(spin,angle,singleparticle,izero,v):
    rad3=sqrt(3.0)
    gamma=angle*pi/180.0
    result=(2*singleparticle-1)*(1.0/izero*a2(angle)-1.0/izero*a1(angle))+2.0*spin*1.0/izero*a1(angle)+v*n2(singleparticle,angle)
    return result
def big1(spin,angle,singleparticle,izero,v):
    minus=-1
    result=(t1(spin,angle,singleparticle,izero)*t2(spin,angle,singleparticle,izero))+8.0*spin*singleparticle*(1.0/izero)*(1.0/izero)*a2(angle)*a3(angle)+t3(spin,angle,singleparticle,izero,v)*t4(spin,angle,singleparticle,izero,v)
    return result*minus
def big2(spin,angle,singleparticle,izero,v):
    result=(t1(spin,angle,singleparticle,izero)*t3(spin,angle,singleparticle,izero,v)-4*spin*singleparticle*(1.0/izero)*(1/izero)*a3(angle)*a3(angle))*(t2(spin,angle,singleparticle,izero)*t4(spin,angle,singleparticle,izero,v)-4*spin*singleparticle*(1.0/izero)*(1/izero)*a2(angle)*a2(angle))
    return result
def omega1(spin,angle,singleparticle,izero,v):
    minus=-1
    result=sqrt(0.5*(minus*big1(spin,angle,singleparticle,izero,v)-sqrt(pow(big1(spin,angle,singleparticle,izero,v),2)-4*big2(spin,angle,singleparticle,izero,v))))
    return result
def omega2(spin,angle,singleparticle,izero,v):
    minus=-1
    result=sqrt(0.5*(minus*big1(spin,angle,singleparticle,izero,v)+sqrt(pow(big1(spin,angle,singleparticle,izero,v),2)-4*big2(spin,angle,singleparticle,izero,v))))
    return result
def hmin(spin,angle,singleparticle,izero,v):
    gamma=angle*pi/180.0
    result=(1.0/izero*a2(angle)+1.0/izero*a3(angle))*0.5*(spin+singleparticle)+pow(spin-singleparticle,2)*1/izero*a1(angle)-v*(2*singleparticle-1)/(singleparticle+1)*sin(gamma+pi/6.0)
    return result
def tsd1(spin,angle,singleparticle,izero,v):
    pi=3.141592654
    result=hmin(spin,angle,singleparticle,izero,v)+0.5*(omega1(spin,angle,singleparticle,izero,v)+omega2(spin,angle,singleparticle,izero,v))
    return result
def fun1(x,a,b):
    result=tsd1(x,17,6.5,a,b)
    return result

params, extras = curve_fit(fun1,mydata.spin,mydata.energy)
print(params)