source: trunk/SRC/ToBeReviewed/TRIANGULATION/definetri_e.pro

;+
;
; @file_comments
;
; @categories
;
; @param INDEX
; It is the tick mark index which starts at 0.
;
; @param NX {in}{required}
; The x dimension array
;
; @param NY {in}{required}
; The y dimension array
;
; @returns
;
; @restrictions
;
; @examples
;
; @history
;
; @version
; $Id$
;
;-
FUNCTION numtri, index, nx, ny
;
  compile_opt idl2, strictarrsubs
;

   y=index/nx
   x=index-y*nx
   numtri = (y NE 0)*(nx-1)*(2*(y-1)+1) + (2-(y EQ (ny-1) OR y EQ (ny-1)))*x

   return, numtri
end
;
;+
;
; @file_comments
; Define a triangulation array like <proidl>TRIANGULATE</proidl> but for a
; E-grid type
;
; @categories
; Make contours with E-grid type
;
; @param NX {in}{required}
; The x dimension array
;
; @param NY {in}{required}
; The y dimension array
;
; @param SINGULAR {in}{optional}
; When singular is undefined all rectangles are cut in using the vertical
; diagonal.
; Singular is a vector which contains the rectangles numbers which are cut in
; using the horizontal diagonal.
; The rectangle number is defined by the index (in a nx*ny vector) of the
; lower-left corner of the rectangle.
;
; @keyword SHIFTED
;
; @returns
; Triangles is a 2d array and is dimensions are 3 and (nx-1)*(ny-1).
; Triangles is defined like in the <proidl>TRIANGULATE</proidl> procedure.
;
; @history
; Sebastien Masson (smlod\@ipsl.jussieu.fr)
;                       June 2001
;
; @version
; $Id$
;
; @todo
; seb: documenter SHIFTED
;
;-
FUNCTION definetri_e, nx, ny, singular, SHIFTED=shifted
;
  compile_opt idl2, strictarrsubs
;
   nx = long(nx)
   ny = long(ny)
   triangles = lonarr(3, 2*(nx-1)*(ny-1))
;
; build the base triangulation with the diamond cut in two triangles
; by the vertical diagonal
;
; first line
   index = lindgen(nx-1)
   trinumber = index
   triangles[0, trinumber] = index
   triangles[1, trinumber] = index+1
   triangles[2, trinumber] = index+(nx+1-shifted)
; last line
   index = (ny-1)*nx+lindgen(nx-1)
   trinumber = numtri(index, nx, ny)
   triangles[0, trinumber] = index
   triangles[1, trinumber] = index+(-nx+((index/nx+1-shifted) MOD 2))
   triangles[2, trinumber] = index+1
; other lines
   if ny GT 2 then begin
      index = lindgen(nx, ny)
      index = index[0:nx-2, 1:ny-2]
      index = index[*]
      oddeven = (index/nx+1-shifted) MOD 2
      trinumber = numtri(index, nx, ny)
      triangles[0, trinumber] = index
      triangles[1, trinumber] = index-nx+oddeven
      triangles[2, trinumber] = index+nx+oddeven
      triangles[0, trinumber+1] = index+nx+oddeven
      triangles[1, trinumber+1] = index-nx+oddeven
      triangles[2, trinumber+1] = index+1
   endif
;
; cut the diamond specified by singular in two triangles
; by the horizontal diagonal
;
   IF keyword_set(singular) then BEGIN
      yindex = singular/nx
      otherline = where(yindex NE 0 AND yindex NE (ny-1))
      if otherline[0] NE -1 then begin
         index = singular[otherline]
         oddeven = (index/nx+1-shifted) MOD 2
         trinumber = numtri(index, nx, ny)
         triangles[0, trinumber] = index
         triangles[1, trinumber] = index-nx+oddeven
         triangles[2, trinumber] = index+1
         triangles[0, trinumber+1] = index
         triangles[1, trinumber+1] = index+1
         triangles[2, trinumber+1] = index+nx+oddeven
      endif

   endif
   return, triangles
end