source: util/src/CalculateSurfaceOfCluster.py.in@ a3ffb44

#!@PYTHON@
#
# This code takes the tesselation in tecplot TRIANGLE format and calculates the surface

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

if (len(sys.argv) < 2):
  print "Usage: "+sys.argv[0]+" <dat>"
  sys.exit(1)

points=[]
def distance( id1, id2 ):
  result = 0
  for i in range(3):
    tmp = points[id1][1+i]-points[id2][1+i]
    result += tmp*tmp
  return math.sqrt(result)

# parse the tecplot dat file (the first three lines are tecplot INFO)
datfile=open(sys.argv[1], "r")
line=datfile.readline()
line=datfile.readline()
line=datfile.readline()
counter=0
for line in datfile:
  if not line.strip():  # end on first empty line, indicating end of points and beginning of triangles
    break
  else:
    fields=line.split()
    points.append([counter, float(fields[0]), float(fields[1]), float(fields[2])])
    counter=counter+1
print "There are %d endpoints of triangles." % (counter)

counter=0
triangles=[]
for line in datfile:
  fields=line.split()
  triangles.append([counter, int(fields[0])-1, int(fields[1])-1, int(fields[2])-1])
  counter=counter+1
print "There are %d triangles." % (counter)

# calculate surface (using Heron's formula)
surface=0
for delta in triangles:
  if delta[1] < len(points) and delta[2] < len(points) and delta[3] < len(points):
    a = distance(delta[1], delta[2])
    b = distance(delta[2], delta[3])
    c = distance(delta[3], delta[1])
    s = (a+b+c)/2.
    print "Triangle %d has sides (%g, %g, %g) and s = %g." % (delta[0], a, b, c, s)
    surface += math.sqrt(s*(s-a)*(s-b)*(s-c))
  else:
    print "Triangle %d consists of points (%d, %d, %d) of which one is unknown." % (delta[0], delta[1], delta[2], delta[3])
print "The calculated surface is %g angstrom^3." % (surface)