"""Support Cartesian points, vectors and operators (dot products,
cross products, etc)."""

import math

class Point:
	"""The Point class represents a point in 3D."""

	def __init__(self, *args):
		if len(args) == 1:
			a = args[0]
			if isinstance(a, Point):
				self.coord = a.coord
			else:
				# assume sequence of some sort
				self.coord = tuple(a)
		elif len(args) == 3:
			self.coord = ( args[0], args[1], args[2] )
		else:
			raise ValueError("Cannot init Point instance")

	def __str__(self):
		return "Point(%8.3f, %8.3f, %8.3f)" % self.coord

	def __add__(self, v):
		"Return the sum of a point and a vector."
		return Point(self.coord[0] + v.coord[0],
					self.coord[1] + v.coord[1],
					self.coord[2] + v.coord[2])

	def __sub__(self, vp):
		"""Return the difference of a point and a vector
		or the vector between two points."""
		if isinstance(vp, Vector):
			return Point(self.coord[0] - vp.coord[0],
					self.coord[1] - vp.coord[1],
					self.coord[2] - vp.coord[2])
		elif isinstance(vp, Point):
			return Vector(self.coord[0] - vp.coord[0],
					self.coord[1] - vp.coord[1],
					self.coord[2] - vp.coord[2])
		else:
			raise ValueError("Cannot subtract point and garbage")

	def __getitem__(self, n):
		"""Return nth component of coordinate."""
		return self.coord[n]

	def sqdistance(self, pt):
		"Return the square of the distance to another point."
		distsq = 0.0
		for i in range(3):
			d = self.coord[i] - pt.coord[i]
			distsq += d * d
		return distsq

	def distance(self, pt):
		"Return distance to another point."
		return math.sqrt(self.sqdistance(pt))

class Vector(Point):
	"""The Vector class represents a direction in 3D.  A Vector
	is not necessarily normalized."""

	def __str__(self):
		return "Vector(%8.3f, %8.3f, %8.3f)" % self.coord

	def __add__(self, v):
		"Return the sum of two vectors."
		return Vector(self.coord[0] + v.coord[0],
					self.coord[1] + v.coord[1],
					self.coord[2] + v.coord[2])

	def __sub__(self, v):
		"Return the difference of two vectors."
		return Vector(self.coord[0] - v.coord[0],
				self.coord[1] - v.coord[1],
				self.coord[2] - v.coord[2])

	def __mul__(self, vs):
		"Return the dot product or scaled vector."
		if isinstance(vs, Vector):
			dot = 0.0
			for i in range(3):
				dot += self.coord[i] * vs.coord[i]
			return dot
		else:
			return Vector(self.coord[0] * vs,
					self.coord[1] * vs,
					self.coord[2] * vs)

	def __rmul__(self, vs):
		"Return the scaled vector."
		return Vector(self.coord[0] * vs,
				self.coord[1] * vs,
				self.coord[2] * vs)

	def __div__(self, scalar):
		"Return the scaled vector."
		return Vector(self.coord[0] / scalar,
				self.coord[1] / scalar,
				self.coord[2] / scalar)

	def normalize(self):
		"Normalize vector (except zero vectors)."
		length = self.length()
		if length != 0.0:
			self.coord = tuple([ c / length for c in self.coord ])

	def sqlength(self):
		"Return the square of the length of the vector."
		lengthsq = 0.0
		for c in self.coord:
			lengthsq += c * c
		return lengthsq

	def length(self):
		"Return the length of the vector."
		return math.sqrt(self.sqlength())

	def cross(self, other):
		"Return the cross product of this vector with other vector"
		u = self.coord
		v = other.coord
		x = u[1] * v[2] - u[2] * v[1]
		y = u[2] * v[0] - u[0] * v[2]
		z = u[0] * v[1] - u[1] * v[0]
		return Vector(x, y, z)

def angle(v0, v1):
	"""Return angle [0..pi] between two vectors"""
	cosa = v0 * v1 / v0.length() / v1.length()
	return math.acos(cosa)

def dihedral(p0, p1, p2, p3):
	"""Return angle [0..2*pi] formed by vertices p0-p1-p2-p3"""
	v01 = p0 - p1
	v32 = p3 - p2
	v12 = p1 - p2
	v0 = v12.cross(v01)
	v3 = v12.cross(v32)
	# The cross product vectors are both normal to the axis
	# vector v12, so the angle between them is the dihedral
	# angle that we are looking for.  However, since "angle"
	# only returns values between 0 and pi, we need to make
	# sure we get the right sign relative to the rotation axis
	a = angle(v0, v3)
	if v0.cross(v3) * v12 > 0:
		a = -a
	return a