;+
; NAME:definetri
;
; PURPOSE:Define a triangulation array like TRIANGULATE. 
;         But in a VERY SIMPLE CASE:
; the points are regulary-gridded on nx*ny array.
; Find a Delaunay triangulation for this set of points is easy:
; Points define (nx-1)*(ny-1) rectangles which we can cut in 2
; triangles. cf. figure above
;
;
;      ny-1*---*---*. . . . . .*---*---*
;          |  /|  /|           |  /|  /|      
;          | / | / |           | / | / |
;          |/  |/  |           |/  |/  |
;      ny-2*---*---*. . . . . .*---*---*   
;          .       .           .       .
;          .       .           .       .
;          .       .           .       .
;         1*---*---*. . . . . .*---*---*
;          |  /|  /|           |  /|  /|
;          | / | / |           | / | / |
;          |/  |/  |           |/  |/  |
;         0*---*---*. . . . . .*---*---*  
;          0   1   2        nx-3  nx-2 nx-1
;
;  You have 2 ways to cut a rectangle:
;      1) the upward diagonal       2) the downward diagonal
;
;          *---*                        *---*
;          |  /|                        |\  |
;          | / |                        | \ |
;          |/  |                        |  \|
;          *---*                        *---*  
;
;
; CATEGORY: to understand how TRIANGULATE and TRIANGULATION work!
;
; CALLING SEQUENCE:triangles=definetri(nx, ny [,downward])
; 
; INPUTS: nx and ny are the array dimensions 
;
; OPTIONAL INPUTS:
;
;        downward: When downward is undefine all rectangles are cut
;        in using the upward diagonal. Downward is a vector which
;        contains the rectangles numbers which are cut in using the
;        downward diagonal.
;        The rectangle number is define by the index (in a nx*ny
;        vector) of the lower-left corner of the rectangle.
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;        triangles is a 2d array and is dimensions are 3 and
;        2*(nx-1)*(ny-1)
;        triangles is define like in the TRIANGULATE procedure.
;        
;
; OPTIONAL OUTPUTS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; EXAMPLE:
;
; triangles=definetri(3,3,[1,3])
; triangles will be a this kind of triangulation:
;
;          *---*---*
;          |\  |  /|
;          | \ | / |
;          |  \|/  |
;          *---*---*
;          |  /|\  |
;          | / | \ |
;          |/  |  \|
;          *---*---*
;
;
; MODIFICATION HISTORY: sebastien Masson (smlod@ipsl.jussieu.fr)
;                       4/3/1999
;
;-
FUNCTION definetri, nx, ny, downward
   nx = long(nx)
   ny = long(ny)
   if n_elements(downward) NE 0 THEN BEGIN 
      if n_elements(downward) GT (nx-1)*(ny-1) then begin
         print, 'downward a trop d''elements par rapport a nx et ny!'
         return,  -1
      endif
      downward = long(downward)
   ENDIF
; we define triangles 
   triangles = lonarr(3, 2*(nx-1)*(ny-1))
;----------------------------------------------------------------------------------
; we cut the rectangles with the upward diagonal
;----------------------------------------------------------------------------------
   if n_elements(downward) NE (nx-1)*(ny-1) then BEGIN ; there is some rectangle to cut.
; we define upward: upward is a vector which contains the rectangles
; numbers which are cut in using the upward diagonal.
; The rectangle number is define by the index (in a nx*ny vector) of
; the lower-left corner of the rectangle.
      upward = bytarr(nx, ny)+1
      upward(*, ny-1) = 0
      upward(nx-1, *) = 0
      if n_elements(downward) NE 0 then upward[downward] = 0
      upward = where(upward EQ 1)
      n1 = n_elements(upward) 
;
; 4 corners indexes of a rectangle number i are
;
;       i+nx  i+nx+1
;          *---*                    
;          |  /|                   
;          | / |                   
;          |/  |                  
;          *---*                
;          i   i+1
;
      trinumber = 2*(upward-upward/nx)
;; we define the right triangles
      triangles[0, trinumber] = upward
      triangles[1, trinumber] = upward+1
      triangles[2, trinumber] = upward+1+nx
; we define the left triangles
      triangles[0, trinumber+1] = upward+1+nx
      triangles[1, trinumber+1] = upward+nx
      triangles[2, trinumber+1] = upward
   ENDIF ELSE n1 = 0
;----------------------------------------------------------------------------------
; we cut the rectangles with the downward diagonal
;----------------------------------------------------------------------------------
   if n_elements(downward) NE 0 then BEGIN
      n2 = n_elements(downward) 
      trinumber = 2*(downward-downward/nx)
; we define the right triangles
      triangles[0, trinumber] = downward+1
      triangles[1, trinumber] = downward+nx+1
      triangles[2, trinumber] = downward+nx
; we define the left triangles
      triangles[0, trinumber+1] = downward+nx
      triangles[1, trinumber+1] = downward
      triangles[2, trinumber+1] = downward+1
   endif

   return, triangles
end