Ignore:
Timestamp:
08/09/06 12:12:54 (18 years ago)
Author:
navarro
Message:

english and nicer header (3a)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/ToBeReviewed/TRIANGULATION/triangule_c.pro

    r134 r150  
    33;------------------------------------------------------------ 
    44;+ 
    5 ; NAME:triangule_c 
    6 ; 
    7 ; PURPOSE:construit le tableau de triangulation. 
    8 ; 
    9 ; L''idee est de 
    10 ; construire une liste de triangles qui relient les points entre 
    11 ; eux. Ceci est fait automatiquement avec la fonction TRIANGULATE. 
    12 ;  ICI: 
    13 ; on tient compte du fait que les points sont disposes sur une grille 
    14 ; (reguliere ou pas, mais pas destructuree, cad que les points sont 
    15 ; ecrits suivant une matrice rectangulaire). Un moyen tres simple de 
    16 ; faire des triangles entre tous les points est alors: 
    17 ; 
    18 ;     pour chaque point (i,j) de la matrice -sauf ceux de la derniere 
    19 ;     ligne et de la derniere colonne- on on appelle le rectangle 
    20 ;     (i,j) le rectangle forme par les 4 points (i,j), (i+1,j), 
    21 ;     (i,j+1), (i+1,j+1). Pour tracer tous les triangles, il suffit de 
    22 ;     tracer les 2 triangles contenus ds les rectangles (i,j) 
    23 ; 
    24 ; au passage on remarque que chaque rectangle (i,j) possede 2 diagonales (si 
    25 ; si faites un dessin c''est vrai), il y a donc 2 choix possibles pour 
    26 ; chaque rectangles qd on veut le couper en 2 triangles... 
     5; 
     6; @file_comments 
     7; Construct the triangulation array. 
     8; 
     9; The idea is: construct a list of triangle which link points between them.  
     10; This is automatically done by the function TRIANGULATE 
     11;  Here: 
     12; we consider the fact that points are disposed on a grid (regular or not,  
     13; but not unstructured, that is to say that points are written following a  
     14; rectangular matrix). A easy way to do triangles between all points is then:  
     15; 
     16;     for each point (i,j) of the matrix -exept those of the last line and of 
     17;     the last column- we call rectangle (i,j) the rectangle made of the four 
     18;     points (i,j), (i+1,j), (i,j+1), (i+1,j+1). To trace all triangle, we just 
     19;     have to trace the 2 triangles contained in rectangles (i,j) 
     20; 
     21; We notice that each rectangle (i,j) have 2 diagonals (it is true... Make a 
     22; drawing to make sure!!), so there are two possible choice for each rectangle 
     23; we want to cut in 2 triangles... 
    2724;  
    28 ; C''est grace a ce choix que l''on va pouvoir tracer les cotes avec 
    29 ; des angles droits. A chaque angle de cote remarquable par 
    30 ; l''existance d''un unique point terre ou d''un unique point mer sur 
    31 ; les 4 cotes d''un rectangle (i,j), il faut couper le rectangle 
    32 ; suivant la diagonale qui qui passe par le point singulier. 
     25; It is thanks to this choice that we will be able to trace coast with right 
     26; angles. At each angle of coast remarkable by the existence of an unique land 
     27; point or of an unique ocean point on one of the four summit of a rectangle (i,j), 
     28; we have to cut the rectangle following the diagonal passing by this point. 
    3329;  
    34 ; CATEGORY:pour faire de beaux graphiques masques 
    35 ; 
    36 ; CALLING SEQUENCE:res=triangule([mask]) 
    37 ; 
    38 ; INPUTS:optionnel:mask c''est le tableau 2d qui sevira a masquer le 
    39 ; champ que l''on tracera apres avec CONTOUR, 
     30; @categories 
     31; graphic 
     32; 
     33; @param MASKENTREE {in}{optional} 
     34; It is a 2d array which will serve to mask the field we will trace after with CONTOUR,  
    4035; ...TRIANGULATION=triangule(mask) 
    41 ; si cet argument n''est pas specifie, la function utilise tmask. 
    42 ; 
    43 ; KEYWORD PARAMETERS: 
    44 ; 
    45 ;       /BASIC: specifie que le masque est sur une grille basice 
    46 ;       (utiliser pour la triangulation ds les coupes verticales et 
    47 ;       des hovmoellers) 
    48 ; 
    49 ;       /KEEP_CONT: to keep the triangulation even on the continents 
    50 ; 
    51 ;       COINMONTE=tableau, pour obtenir le tableau de "coins de terre 
    52 ;       montant" a traiter avec completecointerre.pro ds la variable 
    53 ;       tableau plutot que de la faire passer par la variable globale 
    54 ;       twin_corners_up. 
    55 ; 
    56 ;       COINDESCEND=tableau cf COINMONTE 
    57 ; 
    58 ; OUTPUTS: 
    59 ;       res: tableau 2d (3,nbre de triangles). 
    60 ;    chaque ligne de res represente les indices des points 
    61 ;    constituants les sommets d''un triangle. 
    62 ;    cf. comment on trace les triangles ds dessinetri.pro 
    63 ; 
    64 ; COMMON BLOCKS: 
    65 ;       common.pro different.pro definetri.pro 
    66 ; 
    67 ; SIDE EFFECTS: 
    68 ; 
    69 ; RESTRICTIONS:les donnees dont un veut ensuite faire le contour 
    70 ; doivent etre disposees dans une matrice. Par contre dans la matrice, 
    71 ; la disposition des points peut ne pas etre irreguliere. Si les 
    72 ; donnees sont disposees completement de facon irreguliere, utiliser 
    73 ; TRIANGULE. 
    74 ; 
    75 ; EXAMPLE: 
    76 ; 
    77 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 
     36; If this argument is not specified, the function use tmask 
     37; 
     38; @keyword BASIC 
     39; Specify that the mask is on a basic grid (use the triangulation for vertical cuts and hovmoellers) 
     40; 
     41; @keyword KEEP_CONT 
     42; To keep the triangulation even on the continents 
     43; 
     44; @keyword COINMONTE 
     45; It is an array. To obtain the array of "ascending land corner" to be treated with  
     46; completecointerre.pro in the variable array instead of make it pass by the global  
     47; variable twin_corners_up. 
     48; 
     49; @keyword COINDESCEND 
     50; It is an array. See COINMONTE 
     51; 
     52; @returns 
     53; res: tableau 2d (3,nbre de triangles). 
     54; Each line of res represent indexes of points constituing summits of a triangle.  
     55; See how we trace triangles in definetri.pro 
     56; 
     57; @uses 
     58; common.pro 
     59; different.pro 
     60; definetri.pro 
     61; 
     62; @restrictions 
     63; Datas whose we want to do the contour must be disposed in a matrix.  
     64; On the other hand, in the matrix, the points's arrangement can not be  
     65; irregular. If it is, use TRIANGULE. 
     66; 
     67; @history 
     68; Sebastien Masson (smasson@lodyc.jussieu.fr) 
    7869;                       26/4/1999 
     70; 
     71; @version 
     72; $Id$ 
     73; 
     74; @todo 
     75; seb L.267->268 je ne pense pas que ce soit ce que tu voulais dire mais  
     76; c'est la traduction de ce qu'il y avait écrit. Correction si besoin. 
    7977;- 
    8078;------------------------------------------------------------ 
    … …  
    8381FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, BASIC = basic, KEEP_CONT = keep_cont 
    8482; 
    85   compile_opt idl2, strictarrsubs 
    86 ; 
    87    tempsun = systime(1)         ; pour key_performance 
     83compile_opt idl2, strictarrsubs 
     84; 
     85tempsun = systime(1)            ; For key_performance 
    8886;--------------------------------------------------------- 
    8987@cm_4mesh 
    90   IF NOT keyword_set(key_forgetold) THEN BEGIN 
     88IF NOT keyword_set(key_forgetold) THEN BEGIN 
    9189@updatenew 
    92   ENDIF 
    93 ;------------------------------------------------------------ 
    94 ; le masque est donne ou il faut prendre tmask? 
    95 ;------------------------------------------------------------ 
    96 ; 
    97    msk = maskentree 
    98    taille = size(msk) 
    99    nx = taille[1] 
    100    ny = taille[2] 
    101 ; 
    102    IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular 
    103 ;------------------------------------------------------------ 
    104    if keyword_set(key_periodic)*(nx EQ jpi) $ 
    105     AND NOT keyword_set(basic) then BEGIN  
    106       msk = [msk, msk[0, *]] 
    107       nx = nx+1 
    108    ENDIF 
    109 ;------------------------------------------------------------ 
    110 ; on va trouver la liste des rectangles (i,j) (reperes par leur coin 
    111 ; en bas a gauche) qu''il faut couper suivant une diagonale descendante 
    112 ; on appellera cette liste : pts_downward 
     90ENDIF 
     91;------------------------------------------------------------ 
     92; Is the mask given or do we have to take tmask? 
     93;------------------------------------------------------------ 
     94; 
     95msk = maskentree 
     96taille = size(msk) 
     97nx = taille[1] 
     98ny = taille[2] 
     99; 
     100IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular 
     101;------------------------------------------------------------ 
     102if keyword_set(key_periodic)*(nx EQ jpi) $ 
     103  AND NOT keyword_set(basic) then BEGIN  
     104    msk = [msk, msk[0, *]] 
     105    nx = nx+1 
     106ENDIF 
     107;------------------------------------------------------------ 
     108; We will find the list of rectangles (i,j)(located by their left  
     109; bottom corner) we have to cut folowing a descendant diagonal.  
     110; We will call this list : pts_downward 
    113111;  
    114    pts_downward = 0 
    115  
    116 ; on construit le test qui permet de trouver un tel triangle: 
     112pts_downward = 0 
     113 
     114; We construct the test which allow to find this triangle : 
    117115; 
    118116; 
    … …  
    124122;             msk---------------------shift(msk, -1,  0) 
    125123; 
    126    sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;pts qui entourrent le pt en haut a gauche 
    127    sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;pts qui entourrent le pt en bas a droite 
    128  
    129  
    130    tempdeux = systime(1)        ; pour key_performance =2 
    131 ; pt terre en haut a gauche entoure de pts mer 
    132    liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 ) 
    133    if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
    134 ; pt mer en haut a gauche entoure de pts terre 
    135    liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1) 
    136    if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
    137 ; pt terre en bas a droite entoure de pts mer 
    138    liste = where( (4-sum2)*(1-shift(msk, -1,  0)) EQ 1) 
    139    if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
    140 ; pt mer en bas a droite entoure de pts terre 
    141    liste = where( (1-sum2)*shift(msk, -1,  0) EQ 1) 
    142    if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
    143    undefine, liste 
    144 ; 
    145    IF testvar(var = key_performance) EQ 2 THEN $ 
    146     print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux 
    147 ; 
    148    if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 
    149       tempdeux = systime(1)     ; pour key_performance =2 
    150 ;2 points terre en diagonale montante avec 2 points mer sur la diagonale descendante 
    151       coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $ 
    152                         *(shift(msk, 0, -1)*shift(msk, -1,  0) EQ 1) ) 
    153       if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont] 
    154 ; 
    155       IF testvar(var = key_performance) EQ 2 THEN $ 
    156        print, 'temps triangule: trouver coinmont', systime(1)-tempdeux 
    157       tempdeux = systime(1)     ; pour key_performance =2 
    158 ; 
    159 ;2 points terre en diagonale descendante avec 2 points mer sur la diagonale montante 
    160       coindesc = where( ((1-shift(msk,  0, -1))*(1-shift(msk, -1, 0)) $ 
    161                          *msk*shift(msk, -1, -1) EQ 1) ) 
    162 ; 
    163       IF testvar(var = key_performance) EQ 2 THEN $ 
    164        print, 'temps triangule: trouver coindesc', systime(1)-tempdeux 
    165 ; 
    166     ENDIF 
    167 ; 
    168    if n_elements(pts_downward) EQ 1 then BEGIN  
    169       tempdeux = systime(1)     ; pour key_performance =2 
    170 ; 
    171       triang = definetri(nx, ny) 
    172 ; 
    173       IF testvar(var = key_performance) EQ 2 THEN $ 
    174        print, 'temps triangule: definetri', systime(1)-tempdeux 
    175       coinmont = -1 
    176       coindesc = -1 
    177    ENDIF ELSE BEGIN  
    178       tempdeux = systime(1)     ; pour key_performance =2 
    179       pts_downward = pts_downward[1:n_elements(pts_downward)-1] 
    180       pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))] 
    181 ; aucun rectangle ne peut avoir comme coin en bas a gauche un element 
    182 ; de la derniere colonne ou de la derniere ligne. 
    183 ; il faut donc enlever ces points si ils ont ete selectionnes dans 
    184 ; pts_downward. 
    185       derniere_colonne = (lindgen(ny)+1)*nx-1  
    186       derniere_ligne = lindgen(nx)+(ny-1)*nx  
    187       pts_downward =different(pts_downward,derniere_colonne ) 
    188       pts_downward =different(pts_downward,derniere_ligne ) 
    189       if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 
    190          if coinmont[0] NE -1 then begin 
     124sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1)    ;points which surround the left top point. 
     125sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1)    ;points which surround the right bottom point. 
     126 
     127 
     128tempdeux = systime(1)           ; For key_performance =2 
     129; The left top land point surrounded by ocean points 
     130liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 ) 
     131if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
     132; The left top ocean point surrounded by land points 
     133liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1) 
     134if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
     135; The right bottom land point surrounded by ocean points 
     136liste = where( (4-sum2)*(1-shift(msk, -1,  0)) EQ 1) 
     137if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
     138; The right bottom ocean point surrounded by land points 
     139liste = where( (1-sum2)*shift(msk, -1,  0) EQ 1) 
     140if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 
     141undefine, liste 
     142; 
     143IF testvar(var = key_performance) EQ 2 THEN $ 
     144  print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux 
     145; 
     146if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 
     147    tempdeux = systime(1)       ; For key_performance =2 
     148;2 land points in ascendant diagonal with 2 ocean points in descendant diagonal. 
     149    coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $ 
     150                      *(shift(msk, 0, -1)*shift(msk, -1,  0) EQ 1) ) 
     151    if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont] 
     152; 
     153    IF testvar(var = key_performance) EQ 2 THEN $ 
     154      print, 'temps triangule: trouver coinmont', systime(1)-tempdeux 
     155    tempdeux = systime(1)       ; pour key_performance =2 
     156; 
     157    coindesc = where( ((1-shift(msk,  0, -1))*(1-shift(msk, -1, 0)) $ 
     158                       *msk*shift(msk, -1, -1) EQ 1) ) 
     159; 
     160;2 land points in descendant diagonal with 2 ocean points in ascendant diagonal. 
     161    IF testvar(var = key_performance) EQ 2 THEN $ 
     162      print, 'temps triangule: trouver coindesc', systime(1)-tempdeux 
     163; 
     164ENDIF 
     165; 
     166if n_elements(pts_downward) EQ 1 then BEGIN  
     167    tempdeux = systime(1)       ; For key_performance =2 
     168; 
     169    triang = definetri(nx, ny) 
     170; 
     171    IF testvar(var = key_performance) EQ 2 THEN $ 
     172      print, 'temps triangule: definetri', systime(1)-tempdeux 
     173    coinmont = -1 
     174    coindesc = -1 
     175ENDIF ELSE BEGIN  
     176    tempdeux = systime(1)       ; For key_performance =2 
     177    pts_downward = pts_downward[1:n_elements(pts_downward)-1] 
     178    pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))] 
     179; None rectangle can have an element of the last column or of the  
     180; last line as left bottom corner. 
     181; so we have to remove these points if they has been selected in pts_downward. 
     182    derniere_colonne = (lindgen(ny)+1)*nx-1  
     183    derniere_ligne = lindgen(nx)+(ny-1)*nx  
     184    pts_downward =different(pts_downward,derniere_colonne ) 
     185    pts_downward =different(pts_downward,derniere_ligne ) 
     186    if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 
     187        if coinmont[0] NE -1 then begin 
    191188            coinmont =different(coinmont,derniere_colonne ) 
    192189            coinmont =different(coinmont,derniere_ligne ) 
    193          endif 
    194          if coindesc[0] NE -1 then begin 
     190        endif 
     191        if coindesc[0] NE -1 then begin 
    195192            coindesc =different(coindesc,derniere_colonne ) 
    196193            coindesc =different(coindesc,derniere_ligne ) 
    197          endif 
    198       ENDIF ELSE BEGIN  
    199          coinmont = -1 
    200          coindesc = -1 
    201       ENDELSE  
    202       IF testvar(var = key_performance) EQ 2 THEN $ 
    203        print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux 
    204 ; 
    205       tempdeux = systime(1)     ; pour key_performance =2 
    206       if  pts_downward[0] EQ -1 then triang = definetri(nx, ny) $ 
    207       ELSE triang = definetri(nx, ny, pts_downward) 
    208       IF testvar(var = key_performance) EQ 2 THEN $ 
    209        print, 'temps triangule: definetri', systime(1)-tempdeux 
    210    ENDELSE  
    211 ;------------------------------------------------------------ 
    212 ; on vire les triangles qui ne contiennent que des points terre 
    213 ;------------------------------------------------------------ 
    214 ; 
    215 ;  tres bonne idee qui ne marche pas encore a 200% avec IDL 5.2 
    216 ;  ca devrait aller mieux dans les prochaines versions d''IDL... 
    217 ; 
    218    if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin 
    219       tempdeux = systime(1)     ; pour key_performance =2 
    220 ; on enleve les rectangles qui sont entierement dans la terre 
    221       recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 
    222       IF testvar(var = key_performance) EQ 2 THEN $ 
    223        print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 
    224  
    225 ; en attendant une version qui marche parfaitement, on est contraint 
    226 ; de faire un nouveau tri: 
    227 ; il ne faut pas enlever les rectangles qui n''ont qu''un sommet en 
    228 ; commun. 
     194        endif 
     195    ENDIF ELSE BEGIN  
     196        coinmont = -1 
     197        coindesc = -1 
     198    ENDELSE  
     199    IF testvar(var = key_performance) EQ 2 THEN $ 
     200      print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux 
     201; 
     202    tempdeux = systime(1)       ; For key_performance =2 
     203    if  pts_downward[0] EQ -1 then triang = definetri(nx, ny) $ 
     204    ELSE triang = definetri(nx, ny, pts_downward) 
     205    IF testvar(var = key_performance) EQ 2 THEN $ 
     206      print, 'temps triangule: definetri', systime(1)-tempdeux 
     207ENDELSE  
     208;------------------------------------------------------------ 
     209; We delete land points which only contain land points. 
     210;------------------------------------------------------------ 
     211; 
     212;   
     213if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin 
     214    tempdeux = systime(1)       ; For key_performance =2 
     215; We delete rectangles which are entirely in the land. 
     216    recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 
     217    IF testvar(var = key_performance) EQ 2 THEN $ 
     218      print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 
     219 
     220; We do an other sort : 
     221; We have to do not remove rectangles which only have one common summit. 
    229222; t1 = systime(1) 
    230       indice = intarr(nx, ny) 
    231       trimask = intarr(nx, ny) 
    232       trimask[0:nx-2, 0:ny-2] = 1 
    233       IF recdsterre[0] NE -1 then BEGIN  
    234          tempdeux = systime(1)  ; pour key_performance =2 
    235          indice[recdsterre] = 1 
    236 ;      if NOT keyword_set(basic) then begin 
    237          vire1 = 0 
    238          vire2 = 0 
    239          while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin 
    240 ; vire sont les rectangles qu''il faut retirer de recsterre (en fait 
    241 ; qu''il faut garder bien qu''ils soient entirement dans la terre)   
    242             vire1 = where( (indice*shift(indice, -1, -1) $ 
    243                             *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1) 
    244             if vire1[0] NE -1 THEN BEGIN  
    245                indice[vire1] = 0 
     223    indice = intarr(nx, ny) 
     224    trimask = intarr(nx, ny) 
     225    trimask[0:nx-2, 0:ny-2] = 1 
     226    IF recdsterre[0] NE -1 then BEGIN  
     227        tempdeux = systime(1)   ; For key_performance =2 
     228        indice[recdsterre] = 1 
     229        if NOT keyword_set(basic) then begin 
     230            vire1 = 0 
     231            vire2 = 0 
     232            while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin 
     233; Delete rectangles we have to remove from recsterre (in fact those we have  
     234; to keep although they ar eentirely in the land). 
     235                vire1 = where( (indice*shift(indice, -1, -1) $ 
     236                                *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1) 
     237                if vire1[0] NE -1 THEN BEGIN  
     238                    indice[vire1] = 0 
    246239;               indice[vire1+nx+1] = 0 
    247             endif 
    248              
    249             vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $ 
    250                             *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1) 
    251             if vire2[0] NE -1 THEN BEGIN  
    252                indice[vire2+1] = 0 
     240                endif 
     241                 
     242                vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $ 
     243                                *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1) 
     244                if vire2[0] NE -1 THEN BEGIN  
     245                    indice[vire2+1] = 0 
    253246;               indice[vire2+nx] = 0 
    254             endif 
    255          endwhile 
    256          IF testvar(var = key_performance) EQ 2 THEN $ 
    257           print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux 
    258 ;      endif 
    259          indice[*, ny-1] = 1    ; la deriere colonne te la derniere ligne 
    260          indice[nx-1, *] = 1    ; ne peuvent definir de rectangle. 
    261 ; 
    262          tempdeux = systime(1)  ; pour key_performance =2 
    263          recgarde = where(indice EQ 0) 
    264 ; on recupere les numeros des triangles que l'' on va garder 
    265          trigarde = 2*[recgarde-recgarde/nx] 
    266          trigarde = transpose(temporary(trigarde)) 
    267          trigarde = [trigarde, trigarde+1] 
     247                endif 
     248            endwhile 
     249            IF testvar(var = key_performance) EQ 2 THEN $ 
     250              print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux 
     251        endif 
     252        indice[*, ny-1] = 1     ; The last column and the last line 
     253        indice[nx-1, *] = 1     ; can not define any rectangle. 
     254; 
     255        tempdeux = systime(1)   ; For key_performance =2 
     256        recgarde = where(indice EQ 0) 
     257; We recuperate numbers of triangles we will keep. 
     258        trigarde = 2*[recgarde-recgarde/nx] 
     259        trigarde = transpose(temporary(trigarde)) 
     260        trigarde = [trigarde, trigarde+1] 
    268261;  
    269          triang = triang[*, temporary(trigarde[*])] 
    270          IF testvar(var = key_performance) EQ 2 THEN $ 
     262        triang = triang[*, temporary(trigarde[*])] 
     263        IF testvar(var = key_performance) EQ 2 THEN $ 
    271264          print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux 
    272       endif 
    273    endif 
     265    endif 
     266endif 
    274267; print, 'temps tri triangles', systime(1)-t1  
    275268;------------------------------------------------------------ 
    276 ; quand key_periodic eq 1, triang est une liste d''indice d'un 
    277 ; tableau qui a une colonne de trop. 
    278 ; il faut ramener ca a la matrice initiale en mettant les indivces de 
    279 ; la derniere colonne egaux a ceux de la derniere colonne... 
    280 ;------------------------------------------------------------ 
    281    tempdeux = systime(1)        ; pour key_performance =2 
    282    if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 
    283     AND NOT keyword_set(basic) then BEGIN  
    284       indicey = triang/nx 
    285       indicex = triang-indicey*nx 
    286       nx = nx-1 
    287       liste = where(indicex EQ nx) 
    288       if liste[0] NE -1 then indicex[liste] = 0 
    289       triang = indicex+nx*indicey 
    290       nx = nx+1 
    291       if coinmont[0] NE -1 then begin 
    292          indicey = coinmont/nx 
    293          indicex = coinmont-indicey*nx 
    294          nx = nx-1 
    295          liste = where(indicex EQ nx) 
    296          if liste[0] NE -1 THEN indicex[liste] = 0 
    297          coinmont = indicex+nx*indicey 
    298          nx = nx+1 
    299       endif 
    300       if coindesc[0] NE -1 then begin 
    301          indicey = coindesc/nx 
    302          indicex = coindesc-indicey*nx 
    303          nx = nx-1 
    304          liste = where(indicex EQ nx) 
    305          if liste[0] NE -1 THEN indicex[liste] = 0 
    306          coindesc = indicex+nx*indicey 
    307          nx = nx+1 
    308       endif 
    309    endif 
    310    IF testvar(var = key_performance) EQ 2 THEN $ 
    311     print, 'temps triangule: finitions', systime(1)-tempdeux 
    312  
    313 ;------------------------------------------------------------ 
    314    if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 
    315    if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 
    316 ;------------------------------------------------------------ 
    317   IF NOT keyword_set(key_forgetold) THEN BEGIN 
     269; When key_periodic equal 1, triang is a list of indexes's array which  
     270; have a surplus column. 
     271; We have to put it back to the initial matrix by putting indexes of  
     272; the last column equal to these of the last column... 
     273;------------------------------------------------------------ 
     274tempdeux = systime(1)           ; For key_performance =2 
     275if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 
     276  AND NOT keyword_set(basic) then BEGIN  
     277    indicey = triang/nx 
     278    indicex = triang-indicey*nx 
     279    nx = nx-1 
     280    liste = where(indicex EQ nx) 
     281    if liste[0] NE -1 then indicex[liste] = 0 
     282    triang = indicex+nx*indicey 
     283    nx = nx+1 
     284    if coinmont[0] NE -1 then begin 
     285        indicey = coinmont/nx 
     286        indicex = coinmont-indicey*nx 
     287        nx = nx-1 
     288        liste = where(indicex EQ nx) 
     289        if liste[0] NE -1 THEN indicex[liste] = 0 
     290        coinmont = indicex+nx*indicey 
     291        nx = nx+1 
     292    endif 
     293    if coindesc[0] NE -1 then begin 
     294        indicey = coindesc/nx 
     295        indicex = coindesc-indicey*nx 
     296        nx = nx-1 
     297        liste = where(indicex EQ nx) 
     298        if liste[0] NE -1 THEN indicex[liste] = 0 
     299        coindesc = indicex+nx*indicey 
     300        nx = nx+1 
     301    endif 
     302endif 
     303IF testvar(var = key_performance) EQ 2 THEN $ 
     304  print, 'temps triangule: finitions', systime(1)-tempdeux 
     305 
     306;------------------------------------------------------------ 
     307if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 
     308if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 
     309;------------------------------------------------------------ 
     310IF NOT keyword_set(key_forgetold) THEN BEGIN 
    318311   @updateold 
    319  ENDIF  
    320  
    321    IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun  
    322  
    323    return, triang 
     312ENDIF  
     313 
     314IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun  
     315 
     316return, triang 
    324317 
    325318END  
Note: See TracChangeset for help on using the changeset viewer.