#!@PYTHON@
#
# calculates xray spectra for each time step with coordinates from PDB

import sys, os

exec_prefix="@prefix@"
DATABASE=exec_prefix+"/bin/elements.db"	# contains atomic numbers
XRAY="/media/WORK/tmp/CODICE/XR_Intensity"
PYXPLOT="/home/heber/build/pyxplot/pyxplot"
NAME="ettringite_6Gpa"

# get command line parameters
if (len(sys.argv)) < 5:
	print "Usage: "+sys.argv[0]+" <dir> <prefix> <offsetElement> <offsetXYZ> [start from step]"
	sys.exit(1)
else:
	dir=sys.argv[1]
	prefix=sys.argv[2]
	offsetElement=int(sys.argv[3])
	offsetXYZ=int(sys.argv[4])

# generate hash map: element -> atomic number
atomic={}
input = open(DATABASE,"r")
for line in input:
	if "#" in line:
		continue
	entries=line.split()
	element=entries[1].upper()
	number=int(entries[5])
	print "Putting "+element+" -> "+str(number)+" into hash map."
	atomic[element]=number
input.close()

# parse each PDB
if (len(sys.argv))>=5:
  step=int(sys.argv[5])
else:
  step=0
file=dir+"/"+prefix+("%04d" % (step))+".pdb"
while 1:
	try:
		input=open(file, "r")
	except IOError:
		print file+" is not found, breaking."
		break
	print "Current step is "+str(step)+"."
	# recode pdb to "Symbol atomicnumber X Y Z" structure
	output = open(NAME,"w")
	n=0
	for line in input:
		if "#" in line or "CONECT" in line or "MASTER" in line or "AUTHOR" in line or "COMPND" in line:
			continue
		if "END" in line:
			break
		entries=line.split()
		try:
			element=entries[offsetElement]
		except IndexError:
			print "Could not read element nr."+str(offsetElement)+" in line: "+line+"."
			continue
		number=atomic.get(element.upper(), None)
		if number != None:
			output.write("%s\t%d\t%s\t%s\t%s\n" % (element, int(number), entries[offsetXYZ], entries[offsetXYZ+1], entries[offsetXYZ+2]))
                        n=n+1
		else:
			print "WARNING: atomic number for element "+element+" is not known."
	output.close()
	# call Xray
	outputname=dir+"/xray_intensity"+("%04d" % (step))
	os.system(XRAY+' '+str(n)+' >'+outputname+'.csv')
	# make plotfile and move eps into place
	os.system('/bin/cp -f '+outputname+'.csv C3S_xray_spectrum.csv')
	os.system(PYXPLOT+' xray_spectrum.pyx')
	os.system('/bin/mv C3S_xray_spectrum.eps '+outputname+'.eps')
	# next file
	step+=1
	file=dir+"/"+prefix+("%04d" % (step))+".pdb"

sys.exit(1)