;+
; :Description:
;    Generates a hyperplaneplanar representation of a convex polyhedron.  These planes can then be
;    used to determine if points are inside or outside the polyhedron.
;    This is for a single polyhedron
;
; :Params:
;    vertices = vertices
;    polygons - connectivity of the vertices
;
; :Author: Ronn Kling and Jerry Lefever
;-

pro calcFaceEq, face1X, face1Y, face1Z,a,b,c,d,xc,yc,zc
;calculates the planar equation for each face
;
;Only need three points to define the hyperplane
denom = [[face1X[0:2]], [face1Y[0:2]], [face1Z[0:2]]]
ones = fltarr(3)+1
num1 = [[ones],[face1Y[0:2]], [face1Z[0:2]]]
num2 = [[face1X[0:2]], [ones], [face1Z[0:2]]]
num3 = [[face1X[0:2]], [face1Y[0:2]], [ones]]

ai = determ(num1)/determ(denom)
bi = determ(num2)/determ(denom)
ci = determ(num3)/determ(denom)

denom2 = sqrt(ai^2 + bi^2 + ci^2)

a = ai/denom2
b = bi/denom2
c = ci/denom2
d = -1/denom2

if ( a*xc + b*yc + c*zc + d) lt 0 then begin
  a = -a
  b = -b
  c = -c
  d = -d
endif

return & end

;{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|{{:|

function pointInsideApolyhedron, vertices, polygons, data

;find out how many faces we have
numFaces = n_elements(polygons)/(polygons[0]+1)
numEdges = polygons[0]
;get the center
xc = mean(vertices[0,*])
yc = mean(vertices[1,*])
zc = mean(vertices[2,*])

aArray = fltarr(numFaces)
bArray = fltarr(numFaces)
cArray = fltarr(numFaces)
dArray = fltarr(numFaces)

;get the coefficients for each plane
count = 1L
for i=0,numFaces-1 do begin
  indices = polygons[count:count+numEdges-1]
  currVerts = vertices[*,indices]
  calcFaceEq, reform(currVerts[0,*]), reform(currVerts[1,*]), reform(currVerts[2,*]),a,b,c,d,xc,yc,zc
  aArray[i] = a
  bArray[i] = b
  cArray[i] = c
  dArray[i] = d
  count = count + numEdges + 1
endfor

x = reform(data[0,*])
y = reform(data[1,*])
z = reform(data[2,*])

;go through each point and test where it is
origIndices = lindgen(n_elements(x))
for i=0, numFaces-1 do begin
  Q1 = aArray[i]*x + bArray[i]*y + cArray[i]*z + dArray[i]
  posInd = where( Q1 gt 0, count)
  if count eq 0 then return,-1
  x = x[posInd]
  y = y[posInd]
  z = z[posInd]
  origIndices = origIndices[posInd]
endfor

return, origIndices
end