'''
Angle between domains
'''

from pymol import cmd
import transformations

def centroid_byca(sele):
	center = [0.0, 0.0, 0.0, 0]
	def sumup(x,y,z):
		center[0] += x
		center[1] += y
		center[2] += z
		center[3] += 1
	stored.sumup = sumup
	cmd.iterate_state(-1, 'byca (%s)' % sele, 'stored.sumup(x,y,z)')
	if center[3] == 0:
		return [0.0] * 3
	return [i/center[3] for i in center[:3]]

def domainangle(sele1, sele2, method=cmd.super, quiet=0):
	'''
USAGE

    domainangle sele1, sele2
	'''
	quiet = int(quiet)
	if cmd.is_string(method):
		method = eval(method)

	from math import degrees
	from chempy import cpv

	mobile_tmp = cmd.get_unused_name('tmp')
	cmd.create(mobile_tmp, sele1, zoom=0)

	method(mobile_tmp, sele2)
	mat = cmd.get_object_matrix(mobile_tmp)

	Rv = [mat[i:i+4] for i in [0,4,8,12]]

	cmd.delete(mobile_tmp)

	# calculate the axis & angle
	angle, direction, point = transformations.rotation_from_matrix(Rv)
	center1 = centroid_byca(sele1)
	center2 = centroid_byca(sele2)

	angle = degrees(angle)
	if not quiet:
		print "Angle: %.2f deg, Displacement: %.2f angstrom" % (angle, cpv.distance(center1, center2))
	return angle

cmd.extend('domainangle', domainangle)