# This program will create a 2D NMR Matrix

# import necessary python packages
import math

# tells program where to find Data File
Location = input('Enter Data File Location:')

# opens Data File - "rt" means open for reading text.
Data_File = open(Location,"rt")

# introduces and clears variables
DATA = []
xList = []
yList = []
zList = []
maxList = []

# user inputs parameters
print("Please enter the following parameters.")
Resolution = int(input('Enter Resolution (integer > 0):'))
MPHeight = float(input('Enter Max Peak Height:'))
NRValue = float(input('Enter Noise Reduction Value:'))
Xmin = float(input('Enter minimum F1 Value:'))
Xmax = float(input('Enter maximum F1 Value:'))
Ymin = float(input('Enter minimum F2 Value:'))
Ymax = float(input('Enter maximum F2 Value:'))
GvsL = input('Enter "G" for Gaussian and "L" for Lorentzian:')

# create matrix of loaded data in the format [[x1,y1,z1],[x2,y2,z2],...,[xn,yn,zn]]
with open(Location, "rt") as logFile:
	for line in Data_File:
		line = line[:-1]
		DATA.append(line.split(" "))

# get maximum intensity value for future normalization
n = 0
for n in range(0,len(DATA)):
	if float(Xmin)>float(DATA[n][0]) or float(DATA[n][0])>float(Xmax) or float(Ymin)>float(DATA[n][1]) or float(DATA[n][1])>float(Ymax):
		n+=1
	else:
		maxList.append(float(DATA[n][2]))
maxZ = math.log10(max(maxList))
		
# creates 3 lists while removing any peaks outside of the range
i = 0
for i in range(0,len(DATA)):
	if abs(float(DATA[i][2]))/maxZ < NRValue or float(Xmin) > float(DATA[i][0]) or float(DATA[i][0]) > float(Xmax) or float(Ymin) > float(DATA[i][1]) or float(DATA[i][1]) > float(Ymax):
		i+=1
	else:
		xList.append(float(DATA[i][0]))
		yList.append(float(DATA[i][1]))
		zList.append( abs( MPHeight * math.log10(abs(float( DATA[i][2])))/maxZ ) )

# generate string function for mathematica either using Gaussians or Lorentzians
n = 0
strList = []
print(zList)

for n in range(0,len(xList)):
	Z = str(zList[n])
	W = str(1/float(zList[n]))
	X = str(xList[n])
	Y = str(yList[n])
	# determines in User wanted a Gaussian Function or a Lorentzian Function.
	if GvsL == "G":
		strFunction = "z<"+Z+"*E^(-(5*(x-"+X+"))^2-(5*(y-"+Y+"))^2)"
	else:
		strFunction = "z<"+Z+"*(1+500*"+W+"*(x-"+X+")^2)^-1*(1+500*"+W+"*(y-"+Y+")^2)^-1"
	strList.append(strFunction)
	n+=1

# combines function and prints result
Function = "||".join(strList)
print("RegionPlot3D[z<" + Function + "+0.25,{x,"+str(Xmin)+","+str(Xmax)+"},{y,"+str(Ymin)+","+str(Ymax)+"},{z,0,10},PlotPoints->"+str(Resolution)+"]")

# writes a file on in the directory location
mFunction =  open('Mathematica Function.txt','w')
mFunction.write("RegionPlot3D[" + Function + "+0.25,{x,"+str(Xmin)+","+str(Xmax)+"},{y,"+str(Ymin)+","+str(Ymax)+"},{z,0,10},PlotPoints->"+str(Resolution)+"]")