#!@PYTHON@
#
# calculates pair correlation for a given residuum in a PDB file and the TIP3 molecules

import sys, string, re, math
wrerr = sys.stderr.write
wrout = sys.stdout.write

# check for parameters
if len(sys.argv) < 3:
  wrerr("Usage: "+sys.argv[0]+" <prefix> <offset to xyz> <residuum regexp>\n")
  sys.exit(1)

bin_no=15
offset=int(sys.argv[2])
res=sys.argv[3]
output = open(sys.argv[1]+".min_distance", "w")
output2 = open(sys.argv[1]+".density", "w")
output.write("# timestep\tminimum distance to atom in "+res+" molecule\n")
output2.write("# timestep\t")
for n in range(bin_no):
  output2.write("[%d,%d]\t" % (n,n+1))
output2.write("\n")

timestep=0
TYRMatch=re.compile(res)
TIPMatch=re.compile('TIP3.*WTO')

input = open((sys.argv[1]+".pov1.%04d.pdb" % (timestep)), "r")
while input:
  wrout("Scanning time step %d ... " % (timestep))
  # parse the file
  TYRs=[]
  TIPs=[]
  min=[]
  bins=[]
  for line in input:
    if "#" in line:
      continue
    if "END" in line:
      break
    if "CONECT" in line:
      break
    if TYRMatch.search(line):
      entries=line.split()
      for n in range(3):
        TYRs.append(float(entries[n+offset]))
    if TIPMatch.search(line):
      entries=line.split()
      for n in range(3):
        TIPs.append(float(entries[n+offset]))
  input.close()

  # reset arrays
  center=[0.,0.,0.]
  for j in range(len(TYRs)/3):
    min.append(-1.)
    for n in range(bin_no):
      bins.append(0)
    for n in range(3):
      center[n]+=TYRs[j*3+n]
  for n in range(3):
    center[n]/=len(TYRs)/3
  wrout(" atoms (%d,%d) ... center is at [%f,%f,%f] ... " % (len(TYRs)/3,len(TIPs)/3,center[0], center[1], center[2]))
  # calculate distance
  for i in range(len(TIPs)/3):
    mindistance=-1.
    for j in range(len(TYRs)/3):
      # minimum distance
      distance=0.
      for n in range(3):
        distance+=(TIPs[i*3+n]-TYRs[j*3+n])*(TIPs[i*3+n]-TYRs[j*3+n])
      if (min[j]==-1) or (min[j]*min[j]>distance):  # note down minimum distance
        # radial density
        min[j] = math.sqrt(distance)
        if (mindistance == -1) or (mindistance > min[j]):
          mindistance=min[j]
#   distance=0.
#   for n in range(3):
#     distance+=(TIPs[i*3+n]-center[n])*(TIPs[i*3+n]-center[n])
    for n in range(bin_no):
      if (n*n<mindistance) and ((n+1)*(n+1)>mindistance):
        bins[n]+=1
   # store values
  wrout(" writing ... ")
  output.write("%d\t" % (timestep))
  for n in range(len(TYRs)/3):
    output.write("%lg\t" % (min[n]))
  output.write("\n")
  output2.write("%d\t" % (timestep))
  for n in range(bin_no):
    output2.write("%d\t" % (bins[n]))
  output2.write("\n")
  wrout("done.\n")
  # open next file
  timestep+=1
  input = open((sys.argv[1]+".pov1.%04d.pdb" % (timestep)), "r")