; docformat = 'rst'
;
; NAME: 
;    triangleArea
; PURPOSE:
;    Compute the surface area of a triangle in 3-D space.
;+
; :Description:
;    Compute the surface area of a triangle in 3-D space.
;
; :Categories:
;    Geometry
;
; :Params:
;    v1: in, required, type=fltarr(3)
;        The cartesian coordinates of the first vertex
;    v2: in, required, type=fltarr(3)
;        The cartesian coordinates of the second vertex
;    v3: in, required, type=fltarr(3)
;        The cartesian coordinates of the third vertex
;    v4: out, optional, type=fltarr(3)
;        The cartesian coordinates of the base to v2 (along v1-v3)
;
; :Returns: The surface area of the triangle
;
; :Keywords:
;
; :Examples:
;
; :Uses:
;
; :Author:
;    B.Carry (IMCCE)
;
; :History:
;     Change History::
;        Original Version written in April 2014, B. Carry (IMCCE)
;-
function triangleArea, v1, v2, v3, v4

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  COMPILE_OPT hidden

 ;--I.1-- Input Syntax Verification  ----------------------------------------------------
  if n_params() lt 3 then begin
    message, /IOERROR, 'Syntax -  area = triangleArea(v1, v2, v3)'
    return, -1
  endif


 ;--I.2-- Rewrite Inputs for Readiness  -------------------------------------------------
  x1=v1(0)  &  x2=v2(0)  &  x3=v3(0)
  y1=v1(1)  &  y2=v2(1)  &  y3=v3(1)
  z1=v1(2)  &  z2=v2(2)  &  z3=v3(2)



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Coordinates of the Base-Height Intersection        -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--II.1-- Intersection from Dot-Product and Line Equation  ------------------------------
  dX = x3-x1
  dY = y3-y1
  dZ = z3-z1


 ;--II.2-- General Case  -----------------------------------------------------------------
  if dZ ne 0 then begin

   ;--II.2.1-- Pre-Compute Constants
    C  = z2*dZ
    A1 = x2*dX
    A3 = dX*dX * z1/ dZ
    A4 = x1*dX

    B1 = y2*dY
    B3 = dY*dY * z1/ dZ
    B4 = y1*dY

    D = dZ + dX*dX/dZ + dY*dY/dZ

   ;--II.2.2-- Find z4
    z4 = ( C +A1+A3+A4 + B1+B3+B4 ) / D
    
   ;--II.2.3-- Deduce x4 from line equation
    if dX ne 0 then x4=x1 $
               else x4=x1 + (z4-z1)*dX/dZ

   ;--II.2.4-- Deduce y4 from line equation
    if dY ne 0 then y4=y1 $
               else y4=y1 + (z4-z1)*dY/dZ


 ;--II.3-- Simplification: if Base Aligned with Z  ---------------------------------------
  endif else begin

   ;--II.3.1-- Pre-Compute Constants
    C = x2*dX + y2*dY
    F1= x1 * dY*dY/dX
    F2=-y1 * dY
    D = dX + dY*dY/dX


   ;--II.3.2-- Find z4
    z4 = z1
    
   ;--II.3.3-- Deduce x4 from line equation
    if dX eq 0 then x4=x1 $
               else x4=(C+F1+F2)/D

   ;--II.3.4-- Deduce y4 from line equation
    if dY eq 0 then y4=y1 $
               else y4=y1 + (x4-x1)*dY/dX

  endelse

 ;--II.4-- Coordinates of the Intersection  ----------------------------------------------
  v4=[x4,y4,z4]


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- Surface Area from Simple Triangle Equation         -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--III.1-- Length of Vectors  ------------------------------------------------------------
  base  = sqrt( total( (v3-v1)^2. ) )
  height= sqrt( total( (v4-v2)^2. ) )

 ;--III.2-- Surface Area  -----------------------------------------------------------------
  return, base*height / 2.

end