; plane.pro - Defines a plane surface
;
; ----------------------------------------------------------------------
; CONSTRUCTOR
; ----------------------------------------------------------------------
;Inputs:
;  N - Normal vector to plane. This vector is perpendicular to the plane,
;  and points to the "outside" of the plane.
;  D - May be scalar or vector
;      If scalar: Distance from origin to plane. The closest point on the
;      plane is N*D from the origin
;      If vector: A point on the plane. In this case, the scalar D parameter
;      is calculated such that the plane described passes through this point
function plane::init,N,D
  ; init superclass
  status = self->surface::init()
  self.N=N;
  if n_elements(D) eq 1 then begin
    ;D is really D
    self.D=D;
  end else begin
    ;D is actually r0, a point the plane must pass through.
    self.D=dotp(d,N)
  end
  return, status
end

;function is_inside
;Inputs:
;  R - a vector to test. May be a grid
;Returns:
;  1 if given vector is inside the surface, 0 if not. For a plane, the half
;  of space on the side in the direction of the normal vector is outside,
;  and the half of space on the other side of the plane is inside. Points
;  exactly on the plane are also inside, but don't count on this, due to
;  floating point issues.
function plane::is_inside,r
  return, dotp(r,self.N) le self.D
end

;function intersect
;  Calculates intersection parameter t for this surface and a ray
;  R=R0+V*t
;Inputs:
;  R0 - Origin of ray to test. Must be a single vector
;  V  - Direction of ray to test. May be a grid
;Returns:
;  t parameter of intersection, or negative value or NaN if no intersection.
function plane::intersect,r0,v
  E=dotp(v,self.N)
  result=E*0
  par=where(E eq 0,n_par,complement=npar,n_complement=n_npar)
  if(n_par gt 0) then begin
    if self.D eq dotp(r0,self.N) then result[par]=0 else result[par]=-1
  end
  if(n_npar gt 0) then begin
    result[npar]=(self.D-dotp(r0,self.N))/E[npar]
  end
  return,result
end

;function to_pov
;  Returns a string containing a representation of this surface in POV-Ray code.
;  All vectors in the representation must be converted to left-handed.
function plane::to_pov
  ;The following is correct, it switches handedness to POV (left-handed)
  return,string(self.N[0],self.N[2],self.N[1],self.D,format='(%"plane {  <%25.16e,%25.16e,%25.16e>,%25.16e}")')
end
function plane::to_pov_array
  ;The following is correct, it switches handedness to POV (left-handed)
  return,string(self.N[0],self.N[2],self.N[1],self.D,format='(%"{%25.16e,%25.16e,%25.16e,%25.16e},")')
end

;--------------------------------------------------------------------------
; CLEANUP
; --------------------------------------------------------------------------
; PURPOSE:
;   Lifecycle method to clean itself up. Nothing special in this surface
;   to clean up
;
pro plane::CLEANUP
  self->surface::cleanup
end

; --------------------------------------------------------------------------
; Object Data Structure Definition
; --------------------------------------------------------------------------
; PURPOSE:
;   Define the data structure associated with this object.
pro plane__define
  self = {plane, inherits surface, N:[0d,0d,0d], D:0d}
end