# This program will create a T1 NMR Relaxation Function for Mathematica

# imports necessary python packages
import math

# clears all variables
CSList = []
M0List = []
T1List = []
Continue = "y"

# while user is entering peaks, continue to ask for M0, Chemical Shift, and T1
while Continue == "y":
	ChemicalShift = float(input('Enter Chemical Shift:   '))
	CSList.append(ChemicalShift)
	M0 = input('Enter M0 Value:    ')
	M0List.append(M0)
	T1 = input('Enter T1 Value:    ')
	T1List.append(T1)
	Continue = input('Type y to Continue adding peaks:    ')

# clear's normalized Chemical Shifts list
nCSList = []

# normalize Chemical Shifts to a number between 0 and 1
i = 0
for i in range(0,len(CSList)):
	nCSList.append(CSList[i]/max(CSList))
	i+=1	

xFnctList = []

# make a list of the x function in strings to create lorentzians
i = 0
for i in range(0,len(nCSList)):
	xFnct = "(((90*(x-" + str(nCSList[i]) + "))^2+1)^-1)"
	xFnctList.append(xFnct)
	i += 1
	
yFnctList = []

# make a list of the y function in strings to create relaxation curves
i = 0
for i in range(0,len(nCSList)):
	yFnct = M0List[i] + "*(1-2*E^(-(y*16)/" + T1List[i] + "))"
	yFnctList.append(yFnct)
	i += 1

PeakList = []

# combine the x and y functions to create a 3D function
i = 0
for i in range(0,len(yFnctList)):
	PeakList.append(yFnctList[i] + "*" + xFnctList[i])
	i += 1

# join each of the individual functions
Function = "+".join(PeakList)

# puts the function into a piecewise function to develop a base and 3D model structure
Code = "RegionPlot3D[0.3*Piecewise[{{0,"+Function+"> 0}, {"+Function+","+Function+"<0}}]<z<0.3*Piecewise[{{0,"+Function+"<0},{"+Function+","+Function+">0}}]\[Or]-0.01<z<0.01,{x,-0.2,1.2},{y,-0.02,0.5},{z,-1.2,1.2},PlotPoints -> 200]"

print(Code)