from tkinter import *
from math import *
import scipy
import scipy.optimize
import numpy
import matplotlib
from numpy import arange, sin, pi
matplotlib.use('TkAgg')
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2TkAgg
from matplotlib.figure import Figure


root=Tk()
root.wm_title("MRICalc")

###Magnetisation M0
M=1


#######################################################################
#Pulse selection part

PSframe=LabelFrame(root, text="Pulse sequence selection") 				
PSframe.pack()

def sequenceselection():
   selection = "You selected pulse sequence " + str(var.get())
   label.config(text = selection)


#Buttons of different pulse sequences. Only one choice is possible!
var = IntVar()
R1 = Radiobutton(PSframe, text="Spin-echo", variable=var, value=1,
                  command=sequenceselection)

R1.grid(row=1,column=0)

R2 = Radiobutton(PSframe, text="FLASH", variable=var, value=2,
                  command=sequenceselection)

R2.grid(row=1,column=1)

R3 = Radiobutton(PSframe, text="SSFP-Echo", variable=var, value=3,
                  command=sequenceselection)

R3.grid(row=1,column=2)

R4 = Radiobutton(PSframe, text="TrueFISP", variable=var, value=4,
                  command=sequenceselection)

R4.grid(row=1,column=3)

R5 = Radiobutton(PSframe, text="Dual-echo SSFP", variable=var, value=5,
                  command=sequenceselection)

R5.grid(row=1,column=4)

R6 = Radiobutton(PSframe, text="Dual-echo GRE", variable=var, value=6,
                  command=sequenceselection)

R6.grid(row=1,column=5)

R7 = Radiobutton(PSframe, text="FISP", variable=var, value=7,
                  command=sequenceselection)

R7.grid(row=1,column=6)

label = Label(PSframe)
#label.pack()
label.grid(row=2,column=0, columnspan=3)	

######################################################

# Material selection part

Mframe = LabelFrame(root, text="Material selection")

MLframe = Frame(Mframe)
MLframe.pack(side=LEFT)

MRframe = Frame(Mframe)
MRframe.pack(side=RIGHT)

L = Listbox(MLframe, exportselection = 0)
R = Listbox(MRframe, exportselection = 0)
L.pack(side=LEFT)
R.pack(side=LEFT)

scrollbarL = Scrollbar(MLframe)
scrollbarR = Scrollbar(MRframe)



scrollbarL.pack(side=RIGHT, fill=Y)
scrollbarR.pack(side=RIGHT, fill=Y)


scrollbarL['command'] = L.yview
L['yscrollcommand'] = scrollbarL.set

scrollbarR['command'] = R.yview
R['yscrollcommand'] = scrollbarR.set

#Table of the tissues, order of the tissues is the same as order of the vectors in values
for item in ["Adipose tissue", "Subcutaneous fat", "Marrow fat", "Whole blood (deoxygeneted)", "Whole blood (oxygeneted)", "Cerebrospinal fluid", "Grey matter", "White matter", "Liver", "Mouse liver", "Kidney", "Mouse Kidney", "Muscle", "Front muscle mouse leg ", "Cartilage", "Synovial fluid", "Rat cerebelum", "Rat corpus calosum", "Rat thalamus", "Rat medulla oblongata", "Rat cerebral cortex (right side)", "Rat cerebral cortex (left side)", "Hippocampus (right side)", "Hippocampus (left side)", "Skalice 2", "Skalice 3"]:
    L.insert(END, item)
    R.insert(END, item)

#First value in vector means T1, second value means T2  
values = ((245*10**(-3),70*10**(-3)), (371*10**(-3),133*10**(-3)), (365*10**(-3),133*10**(-3)), (1350*10**(-3),50*10**(-3)), (1350*10**(-3),200*10**(-3)), (4500*10**(-3),2300*10**(-3)), (920*10**(-3),100*10**(-3)), (780*10**(-3),90*10**(-3)), (490*10**(-3),40*10**(-3)), (985*10**(-3),60*10**(-3)), (650*10**(-3),75*10**(-3)), (1259*10**(-3),54*10**(-3)), (1420*10**(-3),31.7*10**(-3)), (1190*10**(-3),39*10**(-3)), (1240*10**(-3),36.9*10**(-3)), (3620*10**(-3),767*10**(-3)), (1343*10**(-3),80*10**(-3)), (1080*10**(-3),128*10**(-3)), (1530*10**(-3),60*10**(-3)), (1560*10**(-3),68*10**(-3)), (1246*10**(-3),76*10**(-3)), (1635*10**(-3),78*10**(-3)), (1249*10**(-3),171*10**(-3)), (1072*10**(-3),254*10**(-3)), (36.59*10**(-3),16.19*10**(-3)), (25.17*10**(-3),15.01*10**(-3)));		#(T1,T2) 


Mframe.pack(side=TOP)



F2 = Frame()
lab = Label(F2)

	

lab.pack()
F2.pack(side=TOP)


#Functions of the first derivations..vector d contains derivations 

## FLASH, Dual-echo GRE
def GREeval(x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	
	d = [M*cos(x[2])*(1-e**(-x[1]/T11))*e**(-x[0]/T21)/(1-e**(-x[1]/T11)*cos(x[2]))-M*sin(x[2])**2*(1-e**(-x[1]/T11))*e**(-x[0]/T21)*e**(-x[1]/T11)/(1-e**(-x[1]/T11)*cos(x[2]))**2-M*cos(x[2])*(1-e**(-x[1]/T12))*e**(-x[0]/T22)/(1-e**(-x[1]/T12)*cos(x[2]))+M*sin(x[2])**2*(1-e**(-x[1]/T12))*e**(-x[0]/T22)*e**(-x[1]/T12)/(1-e**(-x[1]/T12)*cos(x[2]))**2, 
	M*sin(x[2])*e**(-x[1]/T11)*log(e)*e**(-x[0]/T21)/(T11*(1-e**(-x[1]/T11)*cos(x[2])))-M*sin(x[2])*(1-e**(-x[1]/T11))*e**(-x[0]/T21)*e**(-x[1]/T11)*log(e)*cos(x[2])/((1-e**(-x[1]/T11)*cos(x[2]))**2*T11)-M*sin(x[2])*e**(-x[1]/T12)*log(e)*e**(-x[0]/T22)/(T12*(1-e**(-x[1]/T12)*cos(x[2])))+M*sin(x[2])*(1-e**(-x[1]/T12))*e**(-x[0]/T22)*e**(-x[1]/T12)*log(e)*cos(x[2])/((1-e**(-x[1]/T12)*cos(x[2]))**2*T12), 
	-M*sin(x[2])*(1-e**(-x[1]/T11))*e**(-x[0]/T21)*log(e)/((1-e**(-x[1]/T11)*cos(x[2]))*T21)+M*sin(x[2])*(1-e**(-x[1]/T12))*e**(-x[0]/T22)*log(e)/((1-e**(-x[1]/T12)*cos(x[2]))*T22)]
	return d
	
def GRE (x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	
	valueextrem = M*sin(x[2])*(1-e**(-x[1]/T11))*e**(-x[0]/T21)/(1-e**(-x[1]/T11)*cos(x[2]))-M*sin(x[2])*(1-e**(-x[1]/T12))*e**(-x[0]/T22)/(1-e**(-x[1]/T12)*cos(x[2]))
	return valueextrem

##Spin-echo
def Spinechoeval(x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	#[derivation TR, derivation TE]
	d = [M*(e**(-(1/2)*(x[0]-x[1])/T11)*log(e)/T11-e**(-x[0]/T11)*log(e)/T11)*e**(-x[1]/T21)-M*(e**(-(1/2)*(x[0]-x[1])/T12)*log(e)/T12-e**(-x[0]/T12)*log(e)/T12)*e**(-x[1]/T22),
		-M*e**(-(1/2)*(x[0]-x[1])/T11)*log(e)*e**(-x[1]/T21)/T11-M*(1-2*e**(-(1/2)*(x[0]-x[1])/T11)+e**(-x[0]/T11))*e**(-x[1]/T21)*log(e)/T21+M*e**(-(1/2)*(x[0]-x[1])/T12)*log(e)*e**(-x[1]/T22)/T12+M*(1-2*e**(-(1/2)*(x[0]-x[1])/T12)+e**(-x[0]/T12))*e**(-x[1]/T22)*log(e)/T22]
	return d
	
def Spinecho (x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	
	valueextrem = M*(1-2*e**(-(1/2)*(x[0]-x[1])/T11)+e**(-x[0]/T11))*e**(-x[1]/T21)-M*(1-2*e**(-(1/2)*(x[0]-x[1])/T12)+e**(-x[0]/T12))*e**(-x[1]/T22)
	return valueextrem
	
##SSFP-echo
def SSFPechoeval(x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	#[derivation angle, derivation TR]
	d = [M*(1/2+(1/2)*tan((1/2)*x[0])**2)*(1-(1-e**(-x[1]/T11)*cos(x[0]))*(1-(e**(-x[1]/T21))**2)/sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))+M*tan((1/2)*x[0])*(-e**(-x[1]/T11)*sin(x[0])*(1-(e**(-x[1]/T21))**2)/sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2)+(1/2)*(1-e**(-x[1]/T11)*cos(x[0]))*(1-(e**(-x[1]/T21))**2)*((2*(1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0]))))*(e**(-x[1]/T11)*sin(x[0])-(e**(-x[1]/T21))**2*sin(x[0]))+2*(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))*sin(x[0]))/((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2)**(3/2))-M*(1/2+(1/2)*tan((1/2)*x[0])**2)*(1-(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)/sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))-M*tan((1/2)*x[0])*(-sin(x[0])*(1-(e**(-x[1]/T22))**2)/sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2)+(1/2)*(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)*((2*(1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0]))))*(e**(-x[1]/T12)*sin(x[0])-(e**(-x[1]/T22))**2*sin(x[0]))+2*(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))*sin(x[0]))/((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2)**(3/2)),
		M*tan((1/2)*x[0])*(-e**(-x[1]/T11)*log(e)*cos(x[0])*(1-(e**(-x[1]/T21))**2)/(T11*sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))-(2*(1-e**(-x[1]/T11)*cos(x[0])))*(e**(-x[1]/T21))**2*log(e)/(T21*sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))+(1/2)*(1-e**(-x[1]/T11)*cos(x[0]))*(1-(e**(-x[1]/T21))**2)*((2*(1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0]))))*(e**(-x[1]/T11)*log(e)*cos(x[0])/T11+2*(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0]))*log(e)/T21+(e**(-x[1]/T21))**2*e**(-x[1]/T11)*log(e)/T11)+2*(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2*log(e)/T21-2*(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))*(1+cos(x[0]))**2*e**(-x[1]/T11)*log(e)/T11)/((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2)**(3/2))-M*tan((1/2)*x[0])*(e**(-x[1]/T12)*log(e)*(1-(e**(-x[1]/T22))**2)/(T12*sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))-(2*(e**(-x[1]/T12)-cos(x[0])))*(e**(-x[1]/T22))**2*log(e)/(T22*sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))+(1/2)*(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)*((2*(1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0]))))*(e**(-x[1]/T12)*log(e)*cos(x[0])/T12+2*(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0]))*log(e)/T22+(e**(-x[1]/T22))**2*e**(-x[1]/T12)*log(e)/T12)+2*(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2*log(e)/T22-2*(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))*(1+cos(x[0]))**2*e**(-x[1]/T12)*log(e)/T12)/((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2)**(3/2))]
	return d
	
def SSFPecho (x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	
	valueextrem = M*tan((1/2)*x[0])*(1-(1-e**(-x[1]/T11)*cos(x[0]))*(1-(e**(-x[1]/T21))**2)/sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))-M*tan((1/2)*x[0])*(1-(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)/sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))
	return valueextrem

##FISP
def FISPeval(x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	#[derivation angle, derivation TR]
	d = [M*(1/2+(1/2)*tan((1/2)*x[0])**2)*(1-(e**(-x[1]/T11)-cos(x[0]))*(1-(e**(-x[1]/T21))**2)/sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))+M*tan((1/2)*x[0])*(-sin(x[0])*(1-(e**(-x[1]/T21))**2)/sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2)+(1/2)*(e**(-x[1]/T11)-cos(x[0]))*(1-(e**(-x[1]/T21))**2)*((2*(1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0]))))*(e**(-x[1]/T11)*sin(x[0])-(e**(-x[1]/T21))**2*sin(x[0]))+2*(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))*sin(x[0]))/((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2)**(3/2))-M*(1/2+(1/2)*tan((1/2)*x[0])**2)*(1-(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)/sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))-M*tan((1/2)*x[0])*(-sin(x[0])*(1-(e**(-x[1]/T22))**2)/sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2)+(1/2)*(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)*((2*(1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0]))))*(e**(-x[1]/T12)*sin(x[0])-(e**(-x[1]/T22))**2*sin(x[0]))+2*(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))*sin(x[0]))/((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2)**(3/2)),
		M*tan((1/2)*x[0])*(e**(-x[1]/T11)*log(e)*(1-(e**(-x[1]/T21))**2)/(T11*sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))-(2*(e**(-x[1]/T11)-cos(x[0])))*(e**(-x[1]/T21))**2*log(e)/(T21*sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))+(1/2)*(e**(-x[1]/T11)-cos(x[0]))*(1-(e**(-x[1]/T21))**2)*((2*(1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0]))))*(e**(-x[1]/T11)*log(e)*cos(x[0])/T11+2*(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0]))*log(e)/T21+(e**(-x[1]/T21))**2*e**(-x[1]/T11)*log(e)/T11)+2*(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2*log(e)/T21-2*(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))*(1+cos(x[0]))**2*e**(-x[1]/T11)*log(e)/T11)/((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2)**(3/2))-M*tan((1/2)*x[0])*(e**(-x[1]/T12)*log(e)*(1-(e**(-x[1]/T22))**2)/(T12*sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))-(2*(e**(-x[1]/T12)-cos(x[0])))*(e**(-x[1]/T22))**2*log(e)/(T22*sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))+(1/2)*(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)*((2*(1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0]))))*(e**(-x[1]/T12)*log(e)*cos(x[0])/T12+2*(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0]))*log(e)/T22+(e**(-x[1]/T22))**2*e**(-x[1]/T12)*log(e)/T12)+2*(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2*log(e)/T22-2*(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))*(1+cos(x[0]))**2*e**(-x[1]/T12)*log(e)/T12)/((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2)**(3/2))]
	return d
	
def FISP (x, args, args1):
	T11 = float(args[0])
	T21 = float(args[1])
	T12 = float(args1[0])
	T22 = float(args1[1])
	
	valueextrem = M*tan((1/2)*x[0])*(1-(e**(-x[1]/T11)-cos(x[0]))*(1-(e**(-x[1]/T21))**2)/sqrt((1-e**(-x[1]/T11)*cos(x[0])-(e**(-x[1]/T21))**2*(e**(-x[1]/T11)-cos(x[0])))**2-(e**(-x[1]/T21))**2*(1-e**(-x[1]/T11))**2*(1+cos(x[0]))**2))-M*tan((1/2)*x[0])*(1-(e**(-x[1]/T12)-cos(x[0]))*(1-(e**(-x[1]/T22))**2)/sqrt((1-e**(-x[1]/T12)*cos(x[0])-(e**(-x[1]/T22))**2*(e**(-x[1]/T12)-cos(x[0])))**2-(e**(-x[1]/T22))**2*(1-e**(-x[1]/T12))**2*(1+cos(x[0]))**2))
	return valueextrem


################################################
################################################

#Kontrola zda jsou hodnoty  casu a uhlu vetsi nez 0
def valuecheck(x0):
	option=var.get()
	print (x0)
	if (option == 1 and option == 3 and option == 5 and option == 7):
		if (x0[0] <= 0 or x0[1] <= 0 ):
			return 1
		else:
			return 0
	else:
		if (x0[0] <= 0 or x0[1] <= 0 or x0[2] <= 0):
			return 1
		else:
			return 0

def x0value():
	option=var.get()
	if(option == 1 and option == 3 and option == 5 and option == 7):
			q=[0,0]
			
	else:
			q=[0,0,0]
	return q


#################################
#################################	
#Function calculates maximum
def go():
	
	if (len(L.curselection()) and len(R.curselection())):
		idxL = int(L.curselection()[0])
		idxR = int(R.curselection()[0])
		valL = values[idxL]
		valR = values[idxR]
		
		#vectors for fsolve
		vector2 = [1*10**(-3), 1*10**(-3)]
		vector3 = [1, 1*10**(-3), 1*10**(-3)]
		
		
		
		extrem = 0.0
		valext = 0.0
		option=var.get()
		x0=x0value()
		
				
		
		
			
			
		while (valuecheck(x0) == 1):
			if (option == 1):					 
					x0 = scipy.optimize.fsolve(Spinechoeval, vector2, (valL, valR))
					extrem = Spinecho(x0, valL, valR)
					x1 = [x0[0]-x0[0]*0.01,x0[1]-x0[1]*0.01]
					valext = Spinecho(x1, valL, valR)
					vector2 = [vector2[0] + 0.1,vector2[1]+0.01]
					
					#Definition of answer widget
					Answer_repetitiontime_label = Label(Answer, text="Repetition time = ")
					Answer_repetitiontime_value = Label(Answer, text=round(x0[0],3))
					Answer_repetitiontime_unit = Label(Answer, text="s")
						
					Answer_echotime_label = Label(Answer, text="Echo time = ")
					Answer_echotime_value = Label(Answer, text=round(x0[1],3))
					Answer_echotime_unit = Label(Answer, text="s")


					#Position
					Answer_repetitiontime_label.grid(row=0,column=0)
					Answer_repetitiontime_value.grid(row=0,column=1)
					Answer_repetitiontime_unit.grid(row=0,column=2)
						
					Answer_echotime_label.grid(row=1,column=0)
					Answer_echotime_value.grid(row=1,column=1)
					Answer_echotime_unit.grid(row=1,column=2)
										
			if (option == 2):
					x0 = scipy.optimize.fsolve(GREeval, vector3, (valL, valR))
					extrem = GRE(x0, valL, valR)
					x1 = [x0[0]-x0[0]*0.01, x0[1]-x0[1]*0.01, x0[2]-x0[2]*0.01]
					valext = GRE(x1, valL, valR)
					vector3 = [vector3[0] + 0.2,vector3[1]+0.01, vector3[2]+0.01]
					print(x0[1])
					#Definition of answer widget
					Answer_angle_value_label= Label(Answer, text="Angle = ")
					deg=degrees(x0[0]%(2*pi))
					Answer_angle_value = Label(Answer, text=round(deg,1))
					Answer_angle_unit = Label(Answer, text="° ")

					Answer_repetitiontime_label = Label(Answer, text="Repetition time = ")
					Answer_repetitiontime_value = Label(Answer, text=round(x0[1],3))
					Answer_repetitiontime_unit = Label(Answer, text="s")
			
					Answer_echotime_label = Label(Answer, text="Echo time = ")
					Answer_echotime_value = Label(Answer, text=round(x0[2],3))
					Answer_echotime_unit = Label(Answer, text="s")


					#Position
					Answer_angle_value_label.grid(row=0,column=0)
					Answer_angle_value.grid(row=0,column=1)
					Answer_angle_unit.grid(row=0,column=2)

					Answer_repetitiontime_label.grid(row=1,column=0)
					Answer_repetitiontime_value.grid(row=1,column=1)
					Answer_repetitiontime_unit.grid(row=1,column=2)

					Answer_echotime_label.grid(row=2,column=0)
					Answer_echotime_value.grid(row=2,column=1)
					Answer_echotime_unit.grid(row=2,column=2)
					
			if (option == 3):
					x0 = scipy.optimize.fsolve(SSFPechoeval, vector2, (valL, valR))
					extrem = SSFPecho(x00, valL, valR)
					x1 = [x0[0]-x0[0]*0.01,x0[1]-x0[1]*0.01]
					valext = SSFPecho(x1, valL, valR)
					vector2 = [vector2[0] + 0.2,vector2[1]+0.01]
					#Definition of answer widget
					Answer_repetitiontime_label = Label(Answer, text="Repetition time = ")
					Answer_repetitiontime_value = Label(Answer, text=round(x0[0],3))
					Answer_repetitiontime_unit = Label(Answer, text="s")

					Answer_echotime_label = Label(Answer, text="Echo time = ")
					Answer_echotime_value = Label(Answer, text=round(x0[1],3))
					Answer_echotime_unit = Label(Answer, text="s")


					#Position
					Answer_repetitiontime_label.grid(row=0,column=0)
					Answer_repetitiontime_value.grid(row=0,column=1)
					Answer_repetitiontime_unit.grid(row=0,column=2)

					Answer_echotime_label.grid(row=1,column=0)
					Answer_echotime_value.grid(row=1,column=1)
					Answer_echotime_unit.grid(row=1,column=2)
				#if(option == 4):
					#x0 = scipy.optimize.fsolve(func2, vector2, (valL, valR))
					#extrem = GRE(x0, valL, valR)
			if (option == 5):
					x0 = scipy.optimize.fsolve(SSFPechoeval, vector2, (valL, valR))
					extrem = SSFPecho(x0, valL, valR)
					x1 = [x0[0]-x0[0]*0.01,x0[1]-x0[1]*0.01]
					valext = SSFPecho(x1, valL, valR)
					vector2 = [vector2[0] + 0.2,vector2[1]+0.01]
					#Definition of answer widget
					Answer_repetitiontime_label = Label(Answer, text="Repetition time = ")
					Answer_repetitiontime_value = Label(Answer, text=round(x0[0],3))
					Answer_repetitiontime_unit = Label(Answer, text="s")

					Answer_echotime_label = Label(Answer, text="Echo time = ")
					Answer_echotime_value = Label(Answer, text=round(x0[1],3))
					Answer_echotime_unit = Label(Answer, text="s")


					#Position
					Answer_repetitiontime_label.grid(row=0,column=0)
					Answer_repetitiontime_value.grid(row=0,column=1)
					Answer_repetitiontime_unit.grid(row=0,column=2)

					Answer_echotime_label.grid(row=1,column=0)
					Answer_echotime_value.grid(row=1,column=1)
					Answer_echotime_unit.grid(row=1,column=2)
			if (option == 6):
					x0 = scipy.optimize.fsolve(GREeval, vector3, (valL, valR))
					extrem = GRE(x0, valL, valR)
					x1 = [x0[0]-x0[0]*0.01, x0[1]-x0[1]*0.01, x0[2]-x0[2]*0.01]
					valext = GRE(x1, valL, valR)
					vector3 = [vector3[0] + 0.2,vector3[1]+0.01, vector3[2]+0.01]
					#Definition of answer widget
					Answer_angle_value_label= Label(Answer, text="Angle = ")
					deg=degrees(x0[0]%(2*pi))
					Answer_angle_value = Label(Answer, text=round(deg,1))
					Answer_angle_unit = Label(Answer, text="° ")

					Answer_repetitiontime_label = Label(Answer, text="Repetition time = ")
					Answer_repetitiontime_value = Label(Answer, text=round(x0[1],3))
					Answer_repetitiontime_unit = Label(Answer, text="s")
			
					Answer_echotime_label = Label(Answer, text="Echo time = ")
					Answer_echotime_value = Label(Answer, text=round(x0[2],3))
					Answer_echotime_unit = Label(Answer, text="s")


					#Position
					Answer_angle_value_label.grid(row=0,column=0)
					Answer_angle_value.grid(row=0,column=1)
					Answer_angle_unit.grid(row=0,column=2)

					Answer_repetitiontime_label.grid(row=1,column=0)
					Answer_repetitiontime_value.grid(row=1,column=1)
					Answer_repetitiontime_unit.grid(row=1,column=2)

					Answer_echotime_label.grid(row=2,column=0)
					Answer_echotime_value.grid(row=2,column=1)
					Answer_echotime_unit.grid(row=2,column=2)
			if (option == 7):
					x0 = scipy.optimize.fsolve(FISPeval, vector2, (valL, valR))
					extrem = FISP(x00, valL, valR)
					x1 = [x0[0]-x0[0]*0.01,x0[1]-x0[1]*0.01]
					valext = SSFPecho(x1, valL, valR)
					vector2 = [vector2[0] + 0.2,vector2[1]+0.01]
					#Definition of answer widget
					Answer_repetitiontime_label = Label(Answer, text="Repetition time = ")
					Answer_repetitiontime_value = Label(Answer, text=round(x0[0],3))
					Answer_repetitiontime_unit = Label(Answer, text="s")

					Answer_echotime_label = Label(Answer, text="Echo time = ")
					Answer_echotime_value = Label(Answer, text=round(x0[1],3))
					Answer_echotime_unit = Label(Answer, text="s")


					#Position
					Answer_repetitiontime_label.grid(row=0,column=0)
					Answer_repetitiontime_value.grid(row=0,column=1)
					Answer_repetitiontime_unit.grid(row=0,column=2)

					Answer_echotime_label.grid(row=1,column=0)
					Answer_echotime_value.grid(row=1,column=1)
					Answer_echotime_unit.grid(row=1,column=2)
		
		
		
C1 = Button(root, text = "Calculate", command=go,
                 height=3, \
                 width = 12)
C1.pack()

#Definition of answer widget
Answer = LabelFrame(root, text="Answer:", relief=RAISED)
		
Answer.pack()

Answer = LabelFrame(root, text="Answer:", relief=RAISED)
Answer.pack()

############################################
############################################
#GUI Graphs section
Graphs = LabelFrame(root, text="Graphs:", relief=RAISED)
Graphs.pack()

Optionsgraph = LabelFrame(Graphs, text="Options:", relief=RAISED)
Optionsgraph.grid(row=0, column=0)
CheckVar1 = IntVar()
CheckVar2 = IntVar()
CheckVar3 = IntVar()
startinterval = DoubleVar()
stopinterval = DoubleVar()
RTentryval= DoubleVar()
ETentryval= DoubleVar()
FAentryval= DoubleVar()

###Checkbuttons
RTgraph = Checkbutton(Optionsgraph, text = "Repetition time", variable = CheckVar1, \
                 onvalue = 1, offvalue = 0, height=1)
RTgraph.grid(row=0, column=0, ipady=5, sticky=W)

ETgraph = Checkbutton(Optionsgraph, text = "Echo time", variable = CheckVar2, \
                 onvalue = 1, offvalue = 0, height=1)
ETgraph.grid(row=1, column=0, ipady=5, sticky=W)

FAgraph = Checkbutton(Optionsgraph, text = "Flip angle", variable = CheckVar3, \
                 onvalue = 1, offvalue = 0, height=1)
FAgraph.grid(row=2, column=0, ipady=5, sticky=W)

###Entry for prameters of the graph
E1 = Entry(Optionsgraph ,textvariable=startinterval, bd =8, justify=CENTER)
E1.grid(row=3, column=1, ipady=5)

E2 = Entry(Optionsgraph ,textvariable=stopinterval, bd =8, justify=CENTER)
E2.grid(row=3, column=2, ipady=5)

RTentry = Entry(Optionsgraph ,textvariable=RTentryval, bd =8, justify=CENTER)
RTentry.grid(row=0, column=2, ipady=5)

ETentry = Entry(Optionsgraph ,textvariable=ETentryval, bd =8, justify=CENTER)
ETentry.grid(row=1, column=2, ipady=5)

FAentry = Entry(Optionsgraph ,textvariable=FAentryval, bd =8, justify=CENTER)
FAentry.grid(row=2, column=2, ipady=5)


###Labels of fixed values entry
LabelRTentry= Label(Optionsgraph, text="Fixed TR value:")
LabelRTentry.grid(row=0, column=1, ipady=5)

LabelETentry= Label(Optionsgraph, text="Fixed TE value:")
LabelETentry.grid(row=1, column=1, ipady=5)

LabelFAentry= Label(Optionsgraph, text="Fixed flip angle value:")
LabelFAentry.grid(row=2, column=1, ipady=5)

LabelINTentry= Label(Optionsgraph, text="Interval of chosen variable:")
LabelINTentry.grid(row=3, column=0, ipady=5)
#
LabelRTunitsentry= Label(Optionsgraph, text="ms")
LabelRTunitsentry.grid(row=0, column=3, ipady=5)
LabelETunitsentry= Label(Optionsgraph, text="ms")
LabelETunitsentry.grid(row=1, column=3, ipady=5)
LabelFAunitsentry= Label(Optionsgraph, text="°")
LabelFAunitsentry.grid(row=2, column=3, ipady=5)
LabelINTunitsentry= Label(Optionsgraph, text="ms/°")
LabelINTunitsentry.grid(row=3, column=3, ipady=5)


def GREgraph(T11,T12, T21,T22):
	x=[RTentryval.get()/1000,ETentryval.get()/1000,(numpy.pi/180)*FAentryval.get()]
	start=startinterval.get()
	stop=stopinterval.get()
	t = arange(start/1000,stop/1000,0.000001)
	
	if (CheckVar1.get() == 1):
		x[0]=t
	elif (CheckVar2.get() == 1):
		x[1]=t
	elif (CheckVar3.get() == 1):
		x[2]=t
	
	s = abs(M*numpy.sin(x[2])*(1-e**(-x[0]/T11))*e**(-x[1]/T21)/(1-e**(-x[0]/T11)*numpy.cos(x[2]))-M*numpy.sin(x[2])*(1-e**(-x[0]/T12))*e**(-x[1]/T22)/(1-e**(-x[0]/T12)*numpy.cos(x[2])))
	
	
	return s

def Spinechograph(T11,T12, T21,T22):
	x=[RTentryval.get()/1000,ETentryval.get()/1000,(numpy.pi/180)*FAentryval.get()]
	start=startinterval.get()
	stop=stopinterval.get()
	t = arange(start/1000,stop/1000,0.000001)
	
	if (CheckVar1.get() == 1):
		x[0]=t
	elif (CheckVar2.get() == 1):
		x[1]=t
	elif (CheckVar3.get() == 1):
		Displayerror = LabelFrame(Optionsgraph, text="Warning!", relief=RAISED)
		Displayerror.grid(row=4, column=0, columnspan=3)
		Erroranswer = Label(Displayerror, text="Signal is not dependent on flip angle", background="red")
		Erroranswer.grid(row=0, column=0)
	
	s = abs(M*(1-2*e**(-(1/2)*(x[0]-x[1])/T11)+e**(-x[0]/T11))*e**(-x[1]/T21)-M*(1-2*e**(-(1/2)*(x[0]-x[1])/T12)+e**(-x[0]/T12))*e**(-x[1]/T22))
	
	
	return s	
	
def SSFPechograph(T11,T12, T21,T22):
	x=[RTentryval.get()/1000,ETentryval.get()/1000,(numpy.pi/180)*FAentryval.get()]
	start=startinterval.get()
	stop=stopinterval.get()
	t = arange(start/1000,stop/1000,0.000001)
	
	if (CheckVar1.get() == 1):
		x[0]=t
	elif (CheckVar2.get() == 1):
		Displayerror = LabelFrame(Optionsgraph, text="Warning!", relief=RAISED)
		Displayerror.grid(row=4, column=0, columnspan=3)
		Erroranswer = Label(Displayerror, text="Signal is not dependent on echo time", background="red")
		Erroranswer.grid(row=0, column=0)
	elif (CheckVar3.get() == 1):
		x[2]=t
		
	
	s = abs(M*numpy.tan((1/2)*x[2])*(1-(1-e**(-x[0]/T11)*numpy.cos(x[2]))*(1-(e**(-x[0]/T21))**2)/numpy.sqrt((1-e**(-x[0]/T11)*numpy.cos(x[2])-(e**(-x[0]/T21))**2*(e**(-x[0]/T11)-numpy.cos(x[2])))**2-(e**(-x[0]/T21))**2*(1-e**(-x[0]/T11))**2*(1+numpy.cos(x[2]))**2))-M*numpy.tan((1/2)*x[2])*(1-(e**(-x[0]/T12)-numpy.cos(x[2]))*(1-(e**(-x[0]/T22))**2)/numpy.sqrt((1-e**(-x[0]/T12)*numpy.cos(x[2])-(e**(-x[0]/T22))**2*(e**(-x[0]/T12)-numpy.cos(x[2])))**2-(e**(-x[0]/T22))**2*(1-e**(-x[0]/T12))**2*(1+numpy.cos(x[2]))**2)))
	
	return s	

def FISPgraph(T11,T12, T21,T22):
	x=[RTentryval.get()/1000,ETentryval.get()/1000,(numpy.pi/180)*FAentryval.get()]
	start=startinterval.get()
	stop=stopinterval.get()
	t = arange(start/1000,stop/1000,0.000001)
	#t = linspace(start,stop,10)
	if (CheckVar1.get() == 1):
		x[0]=t
	elif (CheckVar2.get() == 1):
		Displayerror = LabelFrame(Optionsgraph, text="Warning!", relief=RAISED)
		Displayerror.grid(row=4, column=0, columnspan=3)
		Erroranswer = Label(Displayerror, text="Signal is not dependent on echo time", background="red")
		Erroranswer.grid(row=0, column=0)
	elif (CheckVar3.get() == 1):
		x[2]=t
		
	
	s = abs(M*numpy.tan((1/2)*x[2])*(1-(e**(-x[0]/T11)-numpy.cos(x[2]))*(1-(e**(-x[0]/T21))**2)/numpy.sqrt((1-e**(-x[0]/T11)*numpy.cos(x[2])-(e**(-x[0]/T21))**2*(e**(-x[0]/T11)-numpy.cos(x[2])))**2-(e**(-x[0]/T21))**2*(1-e**(-x[0]/T11))**2*(1+numpy.cos(x[2]))**2))-M*numpy.tan((1/2)*x[2])*(1-(e**(-x[0]/T12)-numpy.cos(x[2]))*(1-(e**(-x[0]/T22))**2)/numpy.sqrt((1-e**(-x[0]/T12)*numpy.cos(x[2])-(e**(-x[0]/T22))**2*(e**(-x[0]/T12)-numpy.cos(x[2])))**2-(e**(-x[0]/T22))**2*(1-e**(-x[0]/T12))**2*(1+numpy.cos(x[2]))**2)))
	
	return s	

def plotgraph():
	if (len(L.curselection()) and len(R.curselection())):
		idxL = int(L.curselection()[0])
		idxR = int(R.curselection()[0])
		valL = values[idxL]
		valR = values[idxR]
		
	T11 = float(valL[0])
	T21 = float(valL[1])
	T12 = float(valR[0])
	T22 = float(valR[1])
	option=var.get()
	CheckVar11=CheckVar1.get()
	CheckVar12=CheckVar2.get()
	CheckVar13=CheckVar3.get()
	f = Figure(figsize=(6,4), dpi=80)
	a = f.add_subplot(111)
	start=startinterval.get()
	stop=stopinterval.get()
	t = arange(start/1000,stop/1000,0.000001)
	
	
	if (option == 1):					 
			a.plot(Spinechograph(T11,T12, T21,T22))
					
															
	if (option == 2):
			a.plot(GREgraph(T11,T12, T21,T22))
				
										
	if (option == 3):
			a.plot(SSFPechograph(T11,T12, T21,T22))
				
					
				#if(option == 4):
					#x0 = scipy.optimize.fsolve(func2, vector2, (valL, valR))
					#extrem = GRE(x0, valL, valR)
		
	if (option == 5):
			a.plot(SSFPechograph(T11,T12, T21,T22))
				
					
	if (option == 6):
			a.plot(GREgraph(T11,T12, T21,T22))
					
							
	if (option == 7):
			a.plot(FISPgraph(T11,T12, T21,T22))
					
					
	if (CheckVar11 == 1):
		a.set_title('Dependence of the signal on repetition time')
		a.set_xlabel('Repetition time')
	if (CheckVar12 == 1):
		a.set_title('Dependence of the signal on echo time')
		a.set_xlabel('Echo time')
	if (CheckVar13 == 1):
		a.set_title('Dependence of the signal on flip angle')
		a.set_xlabel('Flip angle')
	
	a.set_ylabel('Contrast')

	Displaygraph = LabelFrame(Graphs, text="Display", relief=RAISED)
	Displaygraph.grid(row=0, column=1, rowspan=2)
	canvas = FigureCanvasTkAgg(f, master=Displaygraph)
	canvas.show()
	canvas.get_tk_widget().pack(side=TOP, fill=BOTH, expand=1)

	toolbar = NavigationToolbar2TkAgg( canvas, Displaygraph )
	toolbar.update()
	canvas._tkcanvas.pack(side=TOP, fill=BOTH, expand=1)

 

C2 = Button(Graphs, text = "Plot", command=plotgraph,
                 height=2, \
                 width = 20)
C2.grid(row=1, column=0, sticky=N)



################
################
################
root.mainloop()