source: trunk/SRC/ToBeReviewed/TRIANGULATION/triangule_e.pro @ 231

;+
;
; @file_comments
; Build the triangulation for a E-grid type
;
; @categories
; Graphics
;
; @param MASKENTREE {in}{optional}{type=2d array}
; It is a 2d array which will serve to mask the field we will trace after with CONTOUR,
; ...TRIANGULATION=triangule(mask)
; If this argument is not specified, the function use tmask
;
; @keyword BASIC
; Specify that the mask is on a basic grid (use the triangulation for vertical cuts and hovmoellers)
;
; @keyword COINMONTE {type=array}
; To obtain the array of "ascending land corner" to be treated with
; <pro>completecointerre</pro> in the variable array instead of make it pass 
; by the global variable twin_corners_up.
;
; @keyword COINDESCEND {type=array}
; See COINMONTE
;
; @keyword SHIFTED
;
; @uses
; common.pro
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      june 2001
;
; @version
; $Id$
;
; @todo
; seb L.152->153 je ne pense pas que ce soit ce que tu voulais dire mais
; c'est la traduction de ce qu'il y avait écrit. Correction si besoin.
;-
;
FUNCTION triangule_e, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $
                  , SHIFTED = shifted, BASIC = basic
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;---------------------------------------------------------
   tempsun = systime(1)         ; For key_performance
;------------------------------------------------------------
; Is the mask given or do we have to take tmask?
;------------------------------------------------------------
;
   msk = maskentree
   sizem = size(msk)
   nx = sizem[1]
   ny = sizem[2]
;------------------------------------------------------------
   if keyword_set(key_periodic)*(nx EQ jpi) $
    AND NOT keyword_set(basic) then BEGIN
      msk = [msk, msk[0, *]]
      nx = nx+1
   ENDIF
;
; we will find the diamond that must be cut in two triangle using the
; horizontal diagonal.
;
   index = lindgen(nx, ny)
   index = index[0:nx-2, 1:ny-2]
   if n_elements(shifted) EQ 0 then shifted = 1
   oddeven = (index/nx+1-shifted) MOD 2
   msk1 = msk[index]
   msk2 = msk[index+1]
   sum = msk[index-nx+oddeven]+msk[index+nx+oddeven]
   sum1 = msk2+sum
   sum2 = msk1+sum
;
; horizontal
;
   singularpoint = where((msk1 EQ 0 AND sum1 EQ 3) OR (msk1 EQ 1 AND sum1 EQ 0) $
                         OR (msk2 EQ 0 AND sum2 EQ 3) OR (msk2 EQ 1 AND sum2 EQ 0) $
                         OR (sum EQ 0 AND (msk1+msk2) EQ 2) )
   if singularpoint[0] NE -1 then begin
      horizontal = index[singularpoint]
      triang = definetri_e(nx, ny, horizontal, SHIFTED = shifted)
   ENDIF ELSE triang = definetri_e(nx, ny, SHIFTED = shifted)
;    coinmont = index[where(sum EQ 2 AND (msk1+msk2) EQ 0)]
;    coindesc = index[where(sum EQ 0 AND (msk1+msk2) EQ 2)]
;
; we keep only the triangles which are outside the land
; but for some reasons we will in fact delete the land diamond
;
;    allrecinland = where(sum1+msk1 EQ 0)
;    indexallinland = index[allrecinland]
;    otherrec = (lindgen(nx, ny))[0:nx-2, 1:ny-2]
;    otherrec = different(otherrec, indexallinland)
; ;
;    index = lindgen(nx, ny)
;    index = index[0:nx-3, 2:ny-3]
;    out = inter(index, indexallinland)
;    IF out[0] NE -1 THEN begin
;       out = inter(out+1, indexallinland)
;       IF out[0] NE -1 THEN begin
;          out = out-1
;          oddeven = (out/nx+1-shifted) MOD 2
;          out = inter(out-nx+oddeven, otherrec)
;          IF out[0] NE -1 THEN begin
;             out = inter(out+2*nx, otherrec)
;             IF out[0] NE -1 THEN begin
;                out = out-(nx+((out/nx+shifted) MOD 2))
;             endif
;          endif
;       endif
;    ENDIF
;    help,  out
; ;
;    index = lindgen(nx, ny)
;    index = index[0:nx-3, 2:ny-3]
;    out = inter(index, otherrec)
;    IF out[0] NE -1 THEN begin
;       out = inter(out+1, otherrec)
;       IF out[0] NE -1 THEN begin
;          out = out-1
;          oddeven = (out/nx+1-shifted) MOD 2
;          out = inter(out-nx+oddeven, indexallinland)
;          IF out[0] NE -1 THEN begin
;             out = inter(out+2*nx, indexallinland)
;             IF out[0] NE -1 THEN begin
;                out = out-(nx+((out/nx+shifted) MOD 2))
;             endif
;          endif
;       endif
;    endif
;    help,  out
; ;
;    IF out[0] EQ -1 THEN out = different(indexallinland, out) ELSE out = indexallinland
;    triout = numtri(out, nx, ny)
;    triout = [triout, triout+1]
;    goodtri = lindgen(2*(nx-1)*(ny-1))
;    goodtri = different(goodtri, triout)
;    triang = triang[*, temporary(goodtri)]

; ;
;------------------------------------------------------------
; When key_periodic equal 1, triang is a list of indexes's array which
; have a surplus column.
; We have to put it back to the initial matrix by putting indexes of
; the last column equal to these of the last column...
;------------------------------------------------------------
   tempdeux = systime(1)        ; For key_performance =2
   if keyword_set(key_periodic)*(nx-1 EQ jpi) $
    AND NOT keyword_set(basic) then BEGIN
      indicey = triang/nx
      indicex = triang-indicey*nx
      nx = nx-1
      liste = where(indicex EQ nx)
      if liste[0] NE -1 then indicex[liste] = 0
      triang = indicex+nx*indicey
      nx = nx+1
;       if coinmont[0] NE -1 then begin
;          indicey = coinmont/nx
;          indicex = coinmont-indicey*nx
;          nx = nx-1
;          liste = where(indicex EQ nx)
;          if liste[0] NE -1 THEN indicex[liste] = 0
;          coinmont = indicex+nx*indicey
;          nx = nx+1
;       endif
;       if coindesc[0] NE -1 then begin
;          indicey = coindesc/nx
;          indicex = coindesc-indicey*nx
;          nx = nx-1
;          liste = where(indicex EQ nx)
;          if liste[0] NE -1 THEN indicex[liste] = 0
;          coindesc = indicex+nx*indicey
;          nx = nx+1
;       endif
   endif
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps triangule: finitions', systime(1)-tempdeux

;------------------------------------------------------------
    if arg_present(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont
    if arg_present(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc

   IF NOT keyword_set(key_forgetold) THEN BEGIN
    @updateold
   ENDIF
;------------------------------------------------------------


   IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun

   return, triang

END