import numpy as npy
from GetNodalFunctionsCoeff import  GetNodalFunctionsCoeff

def slope(md,*args):
	"""
	SLOPE - compute the surface slope

	Usage:
		sx,sy,s=slope(md)
		sx,sy,s=slope(md,md.results.TransientSolution(1).Surface)
	"""

	#load some variables (it is much faster if the variables are loaded from md once for all) 
	if m.strcmp(md.mesh.meshtype(),'2Dhorizontal'):
		numberofelements=md.mesh.numberofelements
		numberofnodes=md.mesh.numberofvertices
		index=md.mesh.elements
		x=md.mesh.x ; y=md.mesh.y
	else:
		numberofelements=md.mesh.numberofelements2d
		numberofnodes=md.mesh.numberofvertices2d
		index=md.mesh.elements2d
		x=md.mesh.x2d; y=md.mesh.y2d

	if len(args)==0:
		surf=md.geometry.surface
	elif len(args)==1:
		surf=args[0]
	else:
		raise RuntimeError("slope.py usage error")

	#%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
	alpha,beta=GetNodalFunctionsCoeff(index,x,y)[0:2]

	summation=npy.array([[1],[1],[1]])
	sx=npy.dot(surf[index-1]*alpha,summation).reshape(-1,)
	sy=npy.dot(surf[index-1]*beta,summation).reshape(-1,)

	s=npy.sqrt(sx**2+sy**2)

	if m.strcmp(md.mesh.meshtype(),'3D'):
		sx=project3d(md,'vector',sx,'type','element')
		sy=project3d(md,'vector',sy,'type','element')
		s=npy.sqrt(sx**2+sy**2)

	return (sx,sy,s)
