;+
;
; @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
; completecointerre 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
;
; @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