source: trunk/TRIANGULATION/triangule_e.pro @ 2

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:triangule_e
;
; PURPOSE:buid the triangulation for a E-grid type
;
; CATEGORY:
;
; CALLING SEQUENCE:
; 
; INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      june 2001
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION triangule_e, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $
                  , SHIFTED = shifted, REGULIER = regulier, PERIODIQUE = periodique
   tempsun = systime(1)         ; pour key_performance
@common
   ;;
;------------------------------------------------------------
; le masque est donne ou il faut prendre tmask?
;------------------------------------------------------------
;
   msk = maskentree
   sizem = size(msk)
   nx = sizem[1]
   ny = sizem[2]
;------------------------------------------------------------
   if n_elements(periodique) EQ 0 then periodique = keyword_set(key_periodique)
   if keyword_set(key_periodique) and keyword_set(periodique) $
    AND NOT keyword_set(regulier) 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)]

; ;
;------------------------------------------------------------
; quand key_periodique eq 1, triang est une liste d''indice d'un
; tableau qui a une colonne de trop.
; il faut ramener ca a la matrice initiale en mettant les indivces de
; la derniere colonne egaux a ceux de la derniere colonne...
;------------------------------------------------------------
   tempdeux = systime(1)        ; pour key_performance =2
   if keyword_set(key_periodique) and keyword_set(periodique) $
    AND NOT keyword_set(regulier) 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 cointerremont = coinmont
;    if arg_present(coindescend) THEN coindescend = coindesc ELSE cointerredesc = coindesc
;
; 
;------------------------------------------------------------


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

   return, triang

END