Changeset 29 for trunk


Ignore:
Timestamp:
05/02/06 15:32:01 (18 years ago)
Author:
pinsard
Message:

upgrade of TRIANGULATION according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

Location:
trunk
Files:
1 deleted
13 copied

Legend:

Unmodified
Added
Removed

  • trunk/ToBeReviewed/TRIANGULATION/ciseauxtri.pro

    r27 r29  
    3838;------------------------------------------------------------ 
    3939FUNCTION ciseauxtri, triang, glam, gphi, TOUT = tout, _EXTRA = ex 
    40 @common 
     40;--------------------------------------------------------- 
     41@cm_4mesh 
     42  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     43@updatenew 
     44  ENDIF 
     45;--------------------------------------------------------- 
     46  IF NOT keyword_set(key_periodic) AND NOT keyword_set(key_irregular) $ 
     47    AND NOT (!map.projection LE 7 AND !map.projection NE 0) $ 
     48    AND NOT (!map.projection EQ 14 OR !map.projection EQ 15 $ 
     49             OR !map.projection EQ 18) THEN return, triang 
     50; 
    4151   tempsun = systime(1)         ; pour key_performance 
     52; 
    4253   taille = size(glam) 
    4354   nx = taille[1] 
    … …  
    4556 
    4657   tempdeux = systime(1)        ; pour key_performance =2 
    47    z = convert_coord(glam(*),gphi(*),/data,/to_normal)  
    48    x = z(0, *) 
    49    y = z(1, *) 
     58   z = convert_coord(glam[*],gphi[*],/data,/to_normal)  
     59   x = z[0, *] 
     60   y = z[1, *] 
    5061   tempvar = SIZE(TEMPORARY(z)) ; delete z 
    5162   IF testvar(var = key_performance) EQ 2 THEN $ 
    … …  
    6071    OR !map.projection EQ 14 OR !map.projection EQ 15 OR !map.projection EQ 18 then begin 
    6172      tempdeux = systime(1)     ; pour key_performance =2 
    62  
    63       ind = where(( finite((x*y)[triang[0, *]])$ ; points NaN ds la premiere colonne 
    64                     *finite((x*y)[triang[1, *]]) $ ; points NaN ds la 2eme colonne 
    65                     *finite((x*y)[triang[2, *]]) ) EQ 1) ; points NaN ds la 3eme colonne 
    66  
    67       if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1 
     73; 
     74      test = (x*y)[triang] 
     75      test = finite(temporary(test), /nan) 
     76      test = total(temporary(test), 1) 
     77      ind = where(temporary(test) EQ 0) 
     78; 
     79      if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return,  -1 
     80      trichanged = 1b 
    6881      IF testvar(var = key_performance) EQ 2 THEN $ 
    6982       print, 'temps ciseauxtri: recherche points a NAN', systime(1)-tempdeux 
    7083   endif 
    7184; 
    72  
    73    seuil = 5. 
     85   seuil = 5 < (min([nx, ny])-2) 
    7486; 
    7587; maintenant on vire les triangles dont un des cotes a une taille 
    … …  
    7789; 
    7890 
    79    if keyword_set(key_periodique)  then begin 
    80       tempdeux = systime(1)     ; pour key_performance =2 
    81       ind=where( (abs(x[triang[1,*]]-x[triang[0,*]]) LE (!p.position[2]-!p.position[0])/seuil) $ 
    82                  AND (abs(x[triang[2,*]]-x[triang[1,*]]) LE (!p.position[2]-!p.position[0])/seuil) $ 
    83                  AND (abs(x[triang[0,*]]-x[triang[2,*]]) LE (!p.position[2]-!p.position[0])/seuil) $ 
    84 ; 
    85       AND (abs(y[triang[0,*]]-y[triang[1,*]]) LE (!p.position[3]-!p.position[1])/seuil) $ 
    86        AND (abs(y[triang[1,*]]-y[triang[2,*]]) LE (!p.position[3]-!p.position[1])/seuil) $  
    87        AND (abs(y[triang[2,*]]-y[triang[0,*]]) LE (!p.position[3]-!p.position[1])/seuil) )  
    88  
     91   if keyword_set(key_periodic)  then begin 
     92      tempdeux = systime(1)     ; pour key_performance =2 
     93; 
     94      xtri = x[triang] 
     95      xtri = xtri-shift(xtri, 1, 0) 
     96      testx = abs(temporary(xtri)) GT ((!p.position[2]-!p.position[0])/seuil) 
     97      testx = total(temporary(testx), 1) 
     98; 
     99      ytri = y[triang] 
     100      ytri = ytri-shift(ytri, 1, 0) 
     101      testy = abs(temporary(ytri)) GT ((!p.position[3]-!p.position[1])/seuil) 
     102      testy = total(temporary(testy), 1) 
     103;  
     104      test = temporary(testx)+temporary(testy) 
     105      ind=where(temporary(test) EQ 0)  
     106; 
    89107      IF testvar(var = key_performance) EQ 2 THEN $ 
    90108       print, 'temps ciseauxtri: trouver les triangles trop grands', systime(1)-tempdeux 
    91109; 
    92110      tempdeux = systime(1)     ; pour key_performance =2 
    93       if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1 
     111      if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return,  -1 
     112      trichanged = 1b 
    94113      IF testvar(var = key_performance) EQ 2 THEN $ 
    95114       print, 'temps ciseauxtri: virer les triangles trop grands', systime(1)-tempdeux 
    … …  
    99118; valable qd /TOUT est active 
    100119; 
    101 ;    if keyword_set(key_irregular) then begin 
     120    if keyword_set(key_irregular) then begin 
    102121 
    103 ;       tempdeux = systime(1)     ; pour key_performance =2 
    104 ;       seuil = 0 
    105 ;       seuil = .1*max([!p.position[2]-!p.position[0],!p.position[3]-!p.position[1]]) 
    106 ;       ind=where(x[triang[0,*]] GE !p.position[0]-seuil AND x[triang[0,*]] LE !p.position[2]+seuil $ 
    107 ;                 AND x[triang[1,*]] GE !p.position[0]-seuil AND x[triang[1,*]] LE !p.position[2]+seuil $ 
    108 ;                 AND x[triang[2,*]] GE !p.position[0]-seuil AND x[triang[2,*]] LE !p.position[2]+seuil $ 
    109 ; ; 
    110 ;       AND y[triang[0,*]] GE !p.position[1]-seuil AND y[triang[0,*]] LE !p.position[3]+seuil $ 
    111 ;        AND y[triang[1,*]] GE !p.position[1]-seuil AND y[triang[1,*]] LE !p.position[3]+seuil $ 
    112 ;        AND y[triang[2,*]] GE !p.position[1]-seuil AND y[triang[2,*]] LE !p.position[3]+seuil ) 
    113 ;       IF testvar(var = key_performance) EQ 2 THEN $ 
    114 ;        print, 'temps ciseauxtri: trouver les triangles hors du cadre', systime(1)-tempdeux 
    115 ; ; 
    116 ;       tempdeux = systime(1)     ; pour key_performance =2 
    117 ;       if ind[0] NE -1 then triang = triang[*, ind] ELSE return,  -1 
    118 ;       IF testvar(var = key_performance) EQ 2 THEN $ 
    119 ;        print, 'temps ciseauxtri: virer les triangles hors du cadre', systime(1)-tempdeux 
    120 ;    ENDIF 
     122       tempdeux = systime(1)     ; pour key_performance =2 
     123       xtri = x[triang] 
     124       test1 = xtri GE !p.position[0] 
     125       test2 = xtri LE !p.position[2] 
     126       undefine, xtri 
     127       testx = temporary(test1)*temporary(test2) 
     128       testx = total(temporary(testx), 1) 
     129; 
     130       ytri = y[triang] 
     131       test1 = ytri GE !p.position[1] 
     132       test2 = ytri LE !p.position[3] 
     133       undefine, ytri 
     134       testy = temporary(test1)*temporary(test2) 
     135       testy = total(temporary(testy), 1) 
     136; 
     137       test = temporary(testx)*temporary(testy); 
     138; 
     139       ind=where(temporary(test) NE 0) 
     140; 
     141       if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return,  -1 
     142       trichanged = 1b 
     143       IF testvar(var = key_performance) EQ 2 THEN $ 
     144        print, 'temps ciseauxtri: virer les triangles hors du cadre', systime(1)-tempdeux 
     145    ENDIF 
    121146; 
    122147;        dernier tri 
    123148; 
    124    if n_elements(ind) GT 1 then BEGIN 
     149   if keyword_set(trichanged) then BEGIN 
    125150; 
    126151; il faut verifier que les triangles que l''on garde ne 
    … …  
    149174; ont des indices suivant x qui sont 0 et nx-1. Ils appatienent au 
    150175; rectangles de la colonne nx-1 et non de la colonne 0 
    151       if keyword_set(key_periodique) AND nx EQ jpi then BEGIN 
     176      if keyword_set(key_periodic) AND nx EQ jpi then BEGIN 
    152177      indxmax = indxtriang0 > indxtriang2 
    153178      indxtriang = indxmin + (nx-1)*(indxmin EQ 0 AND indxmax EQ (nx-1)) 
    … …  
    199224; 
    200225   if keyword_set(key_performance) THEN print, 'temps ciseauxtri', systime(1)-tempsun  
     226; 
    201227   return,  triang 
    202228end 

  • trunk/ToBeReviewed/TRIANGULATION/completecointerre.pro

    r27 r29  
    1616; KEYWORD PARAMETERS:  _EXTRA  
    1717; 
     18;        CONT_COLOR: the color of the continent. defaut value is 
     19;        (!d.n_colors - 1) < 255 => white 
     20; 
    1821; OUTPUTS: non 
    1922; 
    … …  
    3235;------------------------------------------------------------ 
    3336;------------------------------------------------------------ 
    34 PRO completecointerre, COINMONTE = coinmonte, COINDESCEND = coindescend, INDICEZOOM = indicezoom $ 
     37PRO draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 
     38@cm_4mesh 
     39; the triangle must not be out of the domain 
     40  IF min(lons, max = maxlon) GE lon1 AND maxlon LE lon2 $ 
     41    AND min(lats, max = maxlat) GE lat1 AND maxlat LE lat2 then BEGIN 
     42; the triangle must not be too big  
     43    z = convert_coord(lons, lats, /data, /to_normal)  
     44    alldist = [(z[0, 2]-z[0, 0])^2 + (z[1, 2]-z[1, 0])^2 $ 
     45               , (z[0, 0]-z[0, 1])^2 + (z[1, 0]-z[1, 1])^2 $ 
     46               , (z[0, 1]-z[0, 2])^2 + (z[1, 1]-z[1, 2])^2] 
     47    IF max(alldist) LT seuil^2 THEN polyfill, lons, lats $ 
     48      , color = cont_color, _extra = ex 
     49    return 
     50  ENDIF  
     51end 
     52;------------------------------------------------------------ 
     53;------------------------------------------------------------ 
     54PRO completecointerre, COINMONTE = coinmonte, COINDESCEND = coindescend $ 
     55                       , CONT_COLOR = cont_color, INDICEZOOM = indicezoom $ 
    3556                       , _extra = ex 
    3657@common 
    3758;------------------------------------------------------------ 
    38    tempsun = systime(1)         ; pour key_performance 
     59;   if NOT keyword_set(coinmonte) then return 
     60;   if NOT keyword_set(coindescend) then return 
     61;   if NOT keyword_set(indicezoom) then return 
     62  tempsun = systime(1)          ; pour key_performance 
    3963;------------------------------------------------------------ 
    4064; definitions des vecteurs coinmont et coindesc 
    4165;------------------------------------------------------------ 
    42    if keyword_set(coinmonte) then coinmont = coinmonte $ 
    43    ELSE coinmont = cointerremont 
    44    if keyword_set(coindescend) then coindesc = coindescend $ 
    45    ELSE coindesc = cointerredesc 
     66  if keyword_set(coinmonte) then coinmont = coinmonte $ 
     67  ELSE coinmont = twin_corners_up 
     68  if keyword_set(coindescend) then coindesc = coindescend $ 
     69  ELSE coindesc = twin_corners_dn 
     70  IF NOT keyword_set(cont_color) THEN cont_color = (!d.n_colors-1) <  255 
    4671;------------------------------------------------------------ 
    4772; definition descoordonnees des points numerotes 1,2,3,4,5,6 cf. les 
    4873; schemas en dessous! 
    4974;------------------------------------------------------------ 
    50    tempdeux = systime(1)        ; pour key_performance =2 
    51    if coinmont[0] NE -1 OR coindesc[0] NE -1 then BEGIN 
    52       if keyword_set(indicezoom) then begin 
    53          long1 = glamv[indicezoom] & lati1 = gphiv[indicezoom] 
    54          long2 = glamu[indicezoom] & lati2 = gphiu[indicezoom] 
    55          long3 = glamf[indicezoom] & lati3 = gphif[indicezoom] 
    56          taille = size(indicezoom) 
    57          nx = taille[1] 
     75  tempdeux = systime(1)         ; pour key_performance =2 
     76  if coinmont[0] NE -1 OR coindesc[0] NE -1 then BEGIN 
     77    if keyword_set(indicezoom) then BEGIN 
     78; if we use key_stide, the t, u, v and f points are no more related to 
     79; the same cell because glamf and gphif has be recomputed to be in the 
     80; middle of two t points 
     81      IF total(key_stride) EQ 3 AND finite(glamv[0]*gphiv[0]) NE 0 THEN BEGIN 
     82        long1 = glamv[indicezoom] & lati1 = gphiv[indicezoom] 
    5883      ENDIF ELSE BEGIN  
    59          long1 = glamv & lati1 = gphiv 
    60          long2 = glamu & lati2 = gphiu 
    61          long3 = glamf & lati3 = gphif 
    62          nx = jpi 
    63       ENDELSE  
    64       long4 = long2  & lati4 = lati2 
    65       long5 = long1 & lati5 = lati1 
    66    endif 
    67    IF testvar(var = key_performance) EQ 2 THEN $ 
     84        long1 = glamt[indicezoom] & lati1 = gphif[indicezoom] 
     85      ENDELSE 
     86      IF total(key_stride) EQ 3 AND finite(glamu[0]*gphiu[0]) NE 0 THEN BEGIN 
     87        long2 = glamu[indicezoom] & lati2 = gphiu[indicezoom] 
     88      ENDIF ELSE BEGIN 
     89        long2 = glamf[indicezoom] & lati2 = gphit[indicezoom] 
     90      ENDELSE 
     91      long3 = glamf[indicezoom] & lati3 = gphif[indicezoom] 
     92    ENDIF ELSE BEGIN  
     93      IF total(key_stride) EQ 3 AND finite(glamv[0]*gphiv[0]) NE 0 THEN BEGIN 
     94        long1 = glamv & lati1 = gphiv 
     95      ENDIF ELSE BEGIN  
     96        long1 = glamt & lati1 = gphif 
     97      ENDELSE 
     98      IF total(key_stride) EQ 3 AND finite(glamu[0]*gphiu[0]) NE 0 THEN BEGIN 
     99        long2 = glamu & lati2 = gphiu 
     100      ENDIF ELSE BEGIN  
     101        long2 = glamf & lati2 = gphit 
     102      ENDELSE 
     103      long3 = glamf & lati3 = gphif 
     104    ENDELSE  
     105; 
     106    nx = (size(long1, /dimensions))[0] 
     107    ny = (size(long1, /dimensions))[1] 
     108    seuil = 5 < (min([nx, ny])-2) 
     109    seuil = min([(!p.position[2]-!p.position[0])/seuil $ 
     110                 , (!p.position[3]-!p.position[1])/seuil]) 
     111; 
     112  ENDIF 
     113; 
     114  IF testvar(var = key_performance) EQ 2 THEN $ 
    68115    print, 'temps completecointerre: positions des points', systime(1)-tempdeux 
    69116; 
    … …  
    75122;                     4     
    76123;     t(i+nx)=1    u(i+nx)       t(i+nx+1)=0 
    77 ;                     |      \ 
    78 ;                     |         \ 
     124;                     |    \ 
     125;                     |        \ 
    79126;         1         3 |            \   5 
    80127;       v(i)---------f(i)------------v(i+1) 
    … …  
    85132; 
    86133; 
    87    if coinmont[0] NE -1 then BEGIN 
    88       tempdeux = systime(1)     ; pour key_performance =2 
    89       for id = 0, n_elements(coinmont)-1 do BEGIN 
    90          i = coinmont[id] 
    91          if long1[i] GE lon1 AND long5[i+1] LE lon2 $ 
    92           AND lati2[i] GE lat1 AND lati4[i+nx] LE lat2 then begin 
    93             polyfill, [long1[i], long2[i], long3[i], long4[i+nx], long5[i+1], long3[i]] $ 
    94              , [lati1[i], lati2[i], lati3[i], lati4[i+nx], lati5[i+1], lati3[i]] $ 
    95              , color = c_cont, _extra = ex 
    96          endif 
    97       endfor 
    98       IF testvar(var = key_performance) EQ 2 THEN $ 
    99        print, 'temps completecointerre: trace de cointerremonte', systime(1)-tempdeux 
    100    endif 
     134  if coinmont[0] NE -1 then BEGIN 
     135    tempdeux = systime(1)       ; pour key_performance =2 
     136    for id = 0, n_elements(coinmont)-1 do BEGIN 
     137      i = coinmont[id] 
     138      ii = i MOD nx 
     139      ij = i/nx 
     140; bottom triangle 
     141      lons = [long1[i], long2[i], long3[i]] 
     142      lats = [lati1[i], lati2[i], lati3[i]] 
     143      draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 
     144; upper triangle 
     145      IF ii NE nx-1 AND ij NE ny-1 THEN BEGIN 
     146        lons = [long3[i], long1[i+1], long2[i+nx]] 
     147        lats = [lati3[i], lati1[i+1], lati2[i+nx]] 
     148        draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 
     149      ENDIF  
     150    ENDFOR  
     151    IF testvar(var = key_performance) EQ 2 THEN $ 
     152      print, 'temps completecointerre: trace de cointerremonte', systime(1)-tempdeux 
     153  ENDIF 
    101154;------------------------------------------------------------ 
    102155; cas coin terre en descendante.: 
    … …  
    115168;      t(i)=0      2 u(i)          t(i+1)=1 
    116169; 
    117    if keyword_set(coindescend) then coindesc = coindescend $ 
    118    ELSE coindesc = cointerredesc 
    119    if coindesc[0] NE -1 then begin 
    120       tempdeux = systime(1)     ; pour key_performance =2 
    121       for id = 0, n_elements(coindesc)-1 do BEGIN 
    122          i = coindesc[id] 
    123          if long1[i] GE lon1 AND long5[i+1] LE lon2 $ 
    124           AND lati2[i] GE lat1 AND lati4[i+nx] LE lat2 then begin 
    125             polyfill, [long1[i], long4[i+nx], long3[i], long2[i], long5[i+1], long3[i]] $ 
    126              , [lati1[i], lati4[i+nx], lati3[i], lati2[i], lati5[i+1], lati3[i]] $ 
    127              , color = c_cont, _extra = ex 
    128          endif 
    129       endfor 
    130       IF testvar(var = key_performance) EQ 2 THEN $ 
    131        print, 'temps completecointerre: trace de cointerredescend', systime(1)-tempdeux 
    132    endif 
     170  if coindesc[0] NE -1 then begin 
     171    tempdeux = systime(1)       ; pour key_performance =2 
     172    for id = 0, n_elements(coindesc)-1 do BEGIN 
     173      i = coindesc[id] 
     174      ii = i MOD nx 
     175      ij = i/nx 
     176      IF ii NE nx-1 AND ij NE ny-1 THEN BEGIN 
     177; left triangle 
     178        lons = [long1[i], long3[i], long2[i+nx]] 
     179        lats = [lati1[i], lati3[i], lati2[i+nx]] 
     180        draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 
     181; right triangle 
     182        lons = [long3[i], long2[i], long1[i+1]] 
     183        lats = [lati3[i], lati2[i], lati1[i+1]] 
     184        draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 
     185      ENDIF 
     186    ENDFOR 
     187    IF testvar(var = key_performance) EQ 2 THEN $ 
     188      print, 'temps completecointerre: trace de cointerredescend', systime(1)-tempdeux 
     189  ENDIF 
    133190 
    134191;------------------------------------------------------------ 
    135    IF keyword_set(key_performance) THEN print, 'temps completecointerre', systime(1)-tempsun  
     192  IF keyword_set(key_performance) THEN print, 'temps completecointerre', systime(1)-tempsun  
    136193;------------------------------------------------------------ 
    137    return 
     194  return 
    138195end 

  • trunk/ToBeReviewed/TRIANGULATION/dessinetri.pro

    r27 r29  
    2222;        triangulate) 
    2323; 
    24 ; KEYWORD PARAMETERS: tous ceux de plots 
     24; KEYWORD PARAMETERS:  
     25; 
     26;        All plots or polyfill keywords. 
     27; 
     28;        WAIT=x. to call wait x second between each triangle draw. 
     29; 
     30;        /ONEBYONE: to draw the triangles one by one 
     31; 
     32;        /FILL: to fill the triangles (using polyfill) instead of plotting them 
     33; 
     34;        CHANGECOLOR=n. to change the color of each traingle. n colors 
     35;        will be used and repeted if necessary. 
     36; 
    2537; 
    2638; OUTPUTS: 
    … …  
    4153;------------------------------------------------------------ 
    4254 
    43 PRO dessinetri, tri, x, y, _extra = ex 
     55PRO dessinetri, tri, x, y, WAIT = wait, ONEBYONE = onebyone, FILL = fill, CHANGECOLOR = changecolor, _extra = ex 
    4456@common 
    4557   tempsun = systime(1)         ; pour key_performance 
    46  
    47    if n_params() EQ 3 then begin 
    48       glam = x 
    49       gphi = y 
     58   a = '' 
     59   if n_params() EQ 3 then BEGIN 
     60     CASE size(x, /n_dimensions)+size(y, /n_dimensions) OF 
     61       2:BEGIN 
     62         nx = n_elements(x)  
     63         ny = n_elements(y)  
     64         glam = x#replicate(1., ny) 
     65         gphi = replicate(1., nx)#y 
     66       END  
     67       4:BEGIN 
     68         glam = x 
     69         gphi = y 
     70       END 
     71       ELSE:BEGIN 
     72         dummy = report('x and y inputs of dessinetri must have the same number of dimensions (1 or 2)') 
     73         return 
     74       END 
     75     ENDCASE 
    5076   ENDIF ELSE BEGIN  
    5177      grille,mask,glam,gphi, tri = tri 
    … …  
    5379      tri = ciseauxtri(tri, glam, gphi) 
    5480   ENDELSE  
    55  
     81; 
     82   IF keyword_set(changecolor) THEN BEGIN 
     83      oldname = !d.name 
     84      if !d.name EQ 'PS' OR !d.name EQ 'Z' then BEGIN  
     85         thisos = strupcase(strmid(!version.os_family, 0, 3)) 
     86         CASE thisOS of 
     87            'MAC': set_plot, thisOS 
     88            'WIN': set_plot, thisOS 
     89            ELSE: set_plot, 'X' 
     90         ENDCASE 
     91         ncolors=(!d.n_colors-1) < 255 
     92         set_plot, oldname 
     93       ENDIF ELSE ncolors=(!d.n_colors-1) < 255 
     94       color = 1+indgen(changecolor)*(ncolors/(changecolor-1)) 
     95     ENDIF ELSE color = 0 
     96     color = color#replicate(1, n_elements(tri)/3/n_elements(color)+1)  
     97; 
    5698   tempdeux = systime(1)         ; pour key_performance =2 
    5799   for i = 0L, n_elements(tri)/3-1 do begin 
    58100      t = [tri[*, i], tri[0, i]] 
    59 ;      plots, glam[t], gphi[t], color = i MOD 255, _extra = ex 
    60 ;      wait, .1 
    61       plots, glam[t], gphi[t], color = 0, _extra = ex 
     101      IF keyword_set(fill) THEN $ 
     102        polyfill, glam[t], gphi[t], color = color[i], _extra = ex $ 
     103      ELSE plots, glam[t], gphi[t], color = color[i], _extra = ex 
     104      IF keyword_set(wait) THEN wait, wait 
     105      IF keyword_set(onebyone) THEN read, a, prompt = 'press a key' 
    62106   ENDFOR 
    63107   IF testvar(var = key_performance) EQ 2 THEN $ 

  • trunk/ToBeReviewed/TRIANGULATION/drawcoast_c.pro

    r27 r29  
    1 PRO drawcoast_c, mask, xf, yf, nx, ny, CONT_THICK = cont_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 
    2 @common 
     1PRO drawcoast_c, mask, xf, yf, nx, ny, COAST_COLOR = coast_color, COAST_THICK = coast_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 
     2;--------------------------------------------------------- 
     3@cm_4mesh 
     4  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     5@updatenew 
     6@updatekwd 
     7  ENDIF 
     8;--------------------------------------------------------- 
    39   tempsun = systime(1)         ; pour key_performance 
    410;--------------------------------------------------------- 
    … …  
    612; on trace les segments verticaux: 
    713; 
    8    if NOT keyword_set(xseuil) then xseuil = 5. 
    9    distanceseuil = (!p.position[2]-!p.position[0])/xseuil 
     14   if NOT keyword_set(yseuil) then yseuil = 5. < (min([nx, ny])-2) 
     15   distanceseuil = (!p.position[3]-!p.position[1])/yseuil 
    1016; liste: liste des points i pourlesquels on va tracer un segment entre 
    1117; le point i,j-1 et i,j 
    1218   tempdeux = systime(1)        ; pour key_performance =2 
    13    liste = where(mask+shift(mask, -1, 0) EQ 1 $ 
    14                  AND (xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2 LE distanceseuil) 
     19   liste = where((mask+shift(mask, -1, 0)) EQ 1 $ 
     20                 AND ((xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2) LE distanceseuil^2) 
    1521   IF liste[0] NE -1 THEN BEGIN 
    1622; on recupere lx et ly qui sont les indices ds un tableau 2d des 
    … …  
    2632          print, 'temps tracecote: determiner liste des points concernes par un trait vertical', systime(1)-tempdeux 
    2733         tempdeux = systime(1)  ; pour key_performance =2 
    28          for pt = 0, n_elements(lx)-1 do BEGIN  
     34         for pt = 0L, n_elements(lx)-1 do BEGIN  
    2935            i = lx[pt] & j = ly[pt] 
    3036            plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $ 
    31              , color=c_cote,thick=cont_thick, /normal, _extra = ex 
     37              , color = coast_color, thick = coast_thick, /normal, _extra = ex 
    3238         endfor 
    3339         IF testvar(var = key_performance) EQ 2 THEN $ 
    … …  
    4349; periodique mais pour le plots  
    4450   tempdeux = systime(1)        ; pour key_performance =2 
    45    if keyword_set(key_periodique) AND nx EQ jpi then begin 
     51   if keyword_set(key_periodic) AND nx EQ jpi then begin 
    4652      mask = [mask, mask[0, *]] 
    4753      xf = [xf, xf[0, *]] 
    … …  
    4955      nx = nx+1 
    5056   ENDIF 
    51    if NOT keyword_set(yseuil) then yseuil = 5. 
    52    distanceseuil = (!p.position[3]-!p.position[1])/yseuil 
    53    liste = where(mask+shift(mask, 0, -1) EQ 1 $ 
    54                  AND (xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2 LE distanceseuil) 
     57   if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2) 
     58   distanceseuil = (!p.position[2]-!p.position[0])/xseuil 
     59   liste = where((mask+shift(mask, 0, -1)) EQ 1 $ 
     60                 AND ((xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2) LE distanceseuil^2) 
    5561   IF liste[0] NE -1 THEN BEGIN 
    5662      ly = liste/nx & lx = temporary(liste)-nx*ly 
    … …  
    6369          print, 'temps tracecote: determiner liste des points concernes par un trait horizontal', systime(1)-tempdeux 
    6470         tempdeux = systime(1)  ; pour key_performance =2 
    65          for pt = 0, n_elements(lx)-1 do BEGIN  
     71         for pt = 0L, n_elements(lx)-1 do BEGIN  
    6672            i = lx[pt] & j = ly[pt] 
    6773            plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $ 
    68              , color=c_cote,thick=cont_thick, /normal, _extra = ex 
     74              , color = coast_color, thick = coast_thick, /normal, _extra = ex 
    6975         endfor 
    7076         IF testvar(var = key_performance) EQ 2 THEN $ 

  • trunk/ToBeReviewed/TRIANGULATION/drawcoast_e.pro

    r27 r29  
    1 PRO drawcoast_e, mask, xf, yf, nx, ny, CONT_THICK = cont_thick, YSEUIL = yseuil, XSEUIL = xseuil, onemore = onemore, _extra = ex 
    2 @common 
     1PRO drawcoast_e, mask, xf, yf, nx, ny, COAST_COLOR = coast_color, COAST_THICK = coast_thick, YSEUIL = yseuil, XSEUIL = xseuil, onemore = onemore, _extra = ex 
     2;--------------------------------------------------------- 
     3@cm_4mesh 
     4  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     5@updatenew 
     6@updatekwd 
     7  ENDIF 
     8;--------------------------------------------------------- 
    39   tempsun = systime(1)         ; pour key_performance 
    410;--------------------------------------------------------- 
    5    if keyword_set(key_periodique) AND nx EQ jpi then begin 
     11   if keyword_set(key_periodic) AND nx EQ jpi then begin 
    612      mask = [mask, mask[0, *]] 
    713      xf = [xf, xf[0, *]] 
    … …  
    1319; 
    1420   if NOT keyword_set(onemore) then onemore = 0 
    15    if NOT keyword_set(xseuil) then xseuil = 5. 
     21   if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2) 
    1622   distanceseuil = (!p.position[2]-!p.position[0])/xseuil 
    1723; liste: liste des points i pourlesquels on va tracer un segment  
    … …  
    2026   indexbis = index-nx+((index/nx+onemore) MOD 2) 
    2127   liste = where(mask[index+1]+mask[indexbis] EQ 1 $ 
    22                  AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil) 
     28                 AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil^2) 
    2329   IF liste[0] NE -1 THEN BEGIN 
    2430      index = index[liste] 
    … …  
    2632      for pt = 0, n_elements(index)-1 do begin 
    2733         plots, [xf[index[pt]], xf[indexbis[pt]]], [yf[index[pt]], yf[indexbis[pt]]] $ 
    28           , color=c_cote,thick=cont_thick, /normal, _extra = ex 
     34           , color = coast_color, thick = coast_thick, /normal, _extra = ex 
    2935      endfor 
    3036   ENDIF 
    … …  
    3238; we plot the borders of the diamond in this sense : / 
    3339; 
    34    if NOT keyword_set(xseuil) then xseuil = 5. 
     40   if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2) 
    3541   distanceseuil = (!p.position[2]-!p.position[0])/xseuil 
    3642; liste: liste des points i pourlesquels on va tracer un segment  
    … …  
    3945   indexbis = index+nx+((index/nx+onemore) MOD 2) 
    4046   liste = where(mask[index+1]+mask[indexbis] EQ 1 $ 
    41                  AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil) 
     47                 AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil^2) 
    4248   IF liste[0] NE -1 THEN BEGIN 
    4349      index = index[liste] 
    … …  
    4551      for pt = 0, n_elements(index)-1 do begin 
    4652         plots, [xf[index[pt]], xf[indexbis[pt]]], [yf[index[pt]], yf[indexbis[pt]]] $ 
    47           , color=c_cote,thick=cont_thick, /normal, _extra = ex 
     53           , color = coast_color, thick = coast_thick, /normal, _extra = ex 
    4854      endfor 
    4955   ENDIF 

  • trunk/ToBeReviewed/TRIANGULATION/section.pro

    r27 r29  
    3232;------------------------------------------------------------ 
    3333 
    34 PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints, BOITE = boite, TYPE = type, DIREC = direc 
    35 @common 
    36 ;  
    37 ;--------------------- 
    38    array = litchamp(field) 
    39 ;--------------------- 
    40 ; definition de boite en fonction de endpoints 
     34PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints $ 
     35             , BOXZOOM = boxzoom, TYPE = type, WDEPTH = wdepth $ 
     36             , DIREC = direc, SHOWBUILD = showbuild, ONLYBOX = onlybox $ 
     37             , _extra = ex 
     38; 
     39;--------------------------------------------------------- 
     40; include common 
     41@cm_4mesh 
     42@cm_4data 
     43@cm_4cal 
     44  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     45@updatenew 
     46@updatekwd 
     47  ENDIF 
     48;-------------------------------------------------------------- 
     49;--------------------- 
     50;--------------------- 
     51; definition de boxzoom en fonction de endpoints 
    4152; puis redefinition du domaine 
    42    boite2d = [min([endpoints[0], endpoints[2]]), max([endpoints[0], endpoints[2]]), min([endpoints[1], endpoints[3]]), max([endpoints[1], endpoints[3]])] 
    43 ; 
    44    minprof = 0 
    45    profdefault = 200 
    46    if n_elements(type) EQ 0 then type = 'nothing' 
    47    Case N_Elements(Boite) OF 
    48       1:boite=[boite2d, minprof, boite[0]] 
    49       2:boite=[boite2d, boite[0],boite[1]] 
    50       Else:$  
    51        if strpos(type, 'z') NE -1 THEN boite=[boite2d, minprof, profdefault] $ 
    52       ELSE boite=[boite2d, prof1, prof2] 
    53    ENDCASE 
    54    if strpos(type, 'z') NE -1 THEN BEGIN 
    55       !y.range = [boite[5], boite[4]] ;on garde les yranges (axe z) avant de changer la boite. 
    56       profmax = boite[5] 
    57       vraiboite = boite 
    58       if vargrid EQ 'W' then gdep = gdepw ELSE gdep = gdept 
    59 ; check with vertical grid limits (nearest level) 
    60       gwork = gdep 
    61 ; check the increse or decrese of the z axis 
    62       IF gwork[1] LE gwork[0] THEN gwork = reverse(gdep, 1) 
    63       niveauprof = where(gwork ge boite[5]) & niveauprof = niveauprof[0] 
    64       if niveauprof NE -1 then boite[5] = gwork[niveauprof]+1 
    65    ENDIF 
    66    domdef, boite, GRILLE=[vargrid], /findalways, _extra = ex 
     53  boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $ 
     54               , min([endpoints[1], endpoints[3]], max = ma13), ma13] 
     55; 
     56  minprof = 0. 
     57  profdefault = 200. 
     58  if n_elements(type) EQ 0 then type = 'nothing' 
     59  Case N_Elements(Boxzoom) OF 
     60    0:localbox = [boxzoom2d, minprof, profdefault] 
     61    1:localbox = [boxzoom2d, minprof, boxzoom[0]] 
     62    2:localbox = [boxzoom2d, boxzoom[0]] 
     63    4:if strpos(type, 'z') NE -1 THEN $ 
     64      localbox = [boxzoom2d, minprof, profdefault] ELSE localbox = boxzoom2d 
     65    5:localbox = [boxzoom2d, minprof, boxzoom[4]] 
     66    6:localbox = [boxzoom2d, boxzoom[4:5]] 
     67    Else:BEGIN 
     68      print, report('Bad definition of the box') 
     69      stop 
     70    END 
     71  ENDCASE 
     72  nelbox = n_elements(localbox) 
     73; 
     74  if keyword_set(wdepth) then grillechoice = [vargrid, 'W'] $ 
     75  ELSE grillechoice = vargrid 
     76  domdef, localbox, GRIDTYPE = grillechoice, /findalways, _extra = ex 
     77  grille, -1, -1, -1, -1, nx, ny 
     78; if less than 10 points where found, we apply domdef over the whole domain  
     79; -> problem... why 10 points as a test value??? 
     80; how can we find a good test value? 
     81  IF nx * ny LE 10 THEN domdef, GRIDTYPE = grillechoice, _extra = ex 
    6782; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef 
    68    lon1 = min([endpoints[0], endpoints[2]]) 
    69    lon2 = max([endpoints[0], endpoints[2]]) 
    70    lat1 = min([endpoints[1], endpoints[3]]) 
    71    lat2 = max([endpoints[1], endpoints[3]]) 
    72 ; on recupere la grille sur la boite 
    73    grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz 
     83  lon1 = min([endpoints[0], endpoints[2]], max = lon2) 
     84  lat1 = min([endpoints[1], endpoints[3]], max = lat2) 
     85; we extend the box along the z axis -> i that way the plot will be drawn 
     86; until its bottom part. 
     87  if strpos(type, 'z') NE -1 THEN BEGIN 
     88;on garde les yranges (axe z) avant de changer la boxzoom. 
     89    !y.range = [localbox[nelbox-1], localbox[nelbox-2]]  
     90    if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN 
     91      firstzw = 0 > (firstzw-1) 
     92      lastzw = (lastzw+1) < (jpk-1) 
     93      nzw = lastzw - firstzw + 1 
     94    ENDIF ELSE BEGIN 
     95      firstzt = 0 > (firstzt-1) 
     96      lastzt = (lastzt+1) < (jpk-1) 
     97      nzt = lastzt - firstzt + 1 
     98    ENDELSE 
     99    IF NOT keyword_set(key_forgetold) THEN BEGIN 
     100     @updateold 
     101   ENDIF  
     102  ENDIF 
     103  !y.range = [localbox[nelbox-1], localbox[nelbox-2]]  
     104; on recupere la grille sur la boxzoom 
     105; 
     106  IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat' 
     107  grille, -1, -1, -1, -1, nx, ny, nz $ 
     108          , firstx, firsty, firstz, lastx, lasty, lastz 
     109; the extend the box in the x and y direction in order to maximise 
     110; the number of triangles used to build the section 
     111  firstx = 0 > (firstx - 1) 
     112  lastx = (lastx + 1) < (jpi -1) 
     113  firsty = 0 > (firsty - 1) 
     114  lasty = (lasty + 1) < (jpj -1) 
     115   
     116  domdef, firstx, lastx, firsty, lasty, firstz, lastz $ 
     117          , /index, gridtype = vargrid 
     118 
     119  IF keyword_set(onlybox) THEN return 
     120; 
     121  grille, mask, glam, gphi, gdep, nx, ny, nz $ 
     122          , firstx, firsty, firstz, lastx, lasty, lastz 
     123 
    74124;--------------------- 
    75125; on definit la triangulation qui va nous permetre de determiner la 
    … …  
    77127; aussi bien que sur la mer. suivant le sens de la section -plutot 
    78128; longitude ou plutot latitude- on definit la facon de trianguler 
    79    if strpos(type, 'x') NE -1 then BEGIN  
    80       downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] 
    81       tri = definetri(nx, ny, (downward)[*])  
    82    ENDIF ELSE tri = definetri(nx, ny) 
     129  if strpos(type, 'x') NE -1 then BEGIN  
     130    downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] 
     131    tri = definetri(nx, ny, (downward)[*])  
     132  ENDIF ELSE tri = definetri(nx, ny) 
     133; If we have an irregular grid that is periodic, then it is possible that 
     134; some of the triangle have a very large size (neighborg points on the 
     135; sphere but far away when doing the projection) and should not be 
     136; taken into account..  
     137  IF keyword_set(key_irregular) AND keyword_set(key_periodic) THEN BEGIN 
     138    glamtri = glam[tri] 
     139    glamtri = abs(glamtri - shift(glamtri, 1, 0)) 
     140    good = temporary(glamtri) LT (10.*max(glam)/nx) 
     141    good = where(total(temporary(good), 1) EQ 3) 
     142    tri = (temporary(tri))[*, temporary(good)] 
     143  ENDIF 
    83144;--------------------- 
    84145; equation de la droite suivant laquelle on fait la section 
    85146;--------------------- 
    86    abc = linearequation(endpoints[0:1], endpoints[2:3]) 
    87    glamtri = glam[tri] 
    88    gphitri = gphi[tri] 
     147  abc = linearequation(endpoints[0:1], endpoints[2:3]) 
     148  glamtri = glam[tri] 
     149  gphitri = gphi[tri] 
    89150; quels sont les points de la triangulation qui sont au dessus et au 
    90151; dessous de la droite ? 
    91    if abc[1] NE 0 THEN $ 
     152  if abc[1] NE 0 THEN $ 
    92153    test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $ 
    93    ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) 
    94  
    95    zero123 = total(test, 1) 
     154  ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) 
     155 
     156  zero123 = total(test, 1) 
    96157; to keep: triangles de la triangulation qui sont a cheval sur la droite. 
    97    tokeep1 = where(zero123 EQ 1) 
    98    tokeep2 = where(temporary(zero123) EQ 2) 
    99    tokeep = [tokeep1, tokeep2] 
    100  
    101    test = test[*, tokeep] 
    102    tri = tri[*, tokeep] 
     158  tokeep1 = where(zero123 EQ 1) 
     159  tokeep2 = where(temporary(zero123) EQ 2) 
     160  tokeep = [tokeep1, tokeep2] 
     161 
     162  test = test[*, tokeep] 
     163  tri = tri[*, tokeep] 
    103164; quel est le sommet du triangle qui est seul d''un cote de la droite? 
    104    single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) 
    105    single1 = single1-(single1/3)*3 
    106    single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0) 
    107    single2 = single2-(single2/3)*3 
    108  
    109    undefine, tokeep 
    110    undefine, tokeep1 
    111    undefine, tokeep2 
    112    undefine, test 
    113  
    114    single = [temporary(single1), temporary(single2)] 
     165  single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) 
     166  single1 = single1-(single1/3)*3 
     167  single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0) 
     168  single2 = single2-(single2/3)*3 
     169 
     170  undefine, tokeep 
     171  undefine, tokeep1 
     172  undefine, tokeep2 
     173  undefine, test 
     174 
     175  single = [temporary(single1), temporary(single2)] 
    115176; points1 le point du triangles qui est seul d''un cote de la droite. 
    116177; point2 l''autre point du triangle de l''autre cote de la droite 
    117    point1 = [single, single] 
    118    point2 = [single EQ 0, 1 + (single LE 1)] 
    119  
    120    undefine,  single 
    121  
    122    ntri = (size(tri))[2] 
    123    index = [lindgen(ntri), lindgen(ntri)] 
    124  
    125    points1 = tri[point1, index] 
    126    points2 = tri[point2, temporary(index)] 
     178  point1 = [single, single] 
     179  point2 = [single EQ 0, 1 + (single LE 1)] 
     180 
     181  undefine,  single 
     182 
     183  ntri = (size(tri))[2] 
     184  index = [lindgen(ntri), lindgen(ntri)] 
     185 
     186  points1 = tri[point1, index] 
     187  points2 = tri[point2, temporary(index)] 
    127188; points : complexe contenant les couples de points de part et 
    128189; d''autre de la droite. Ils faut supprimer les doublons. 
    129    points = dcomplex(points1, points2) 
    130    points = points[uniq(points, sort(points))] 
    131    symetrique = dcomplex(imaginary(points), double(points)) 
    132    points = points[where(points-shift(temporary(symetrique), 1) NE 0)] 
     190  points = dcomplex(points1, points2) 
     191  points = points[uniq(points, sort(points))] 
     192  symetrique = dcomplex(imaginary(points), double(points)) 
     193  points = points[where(points-shift(temporary(symetrique), 1) NE 0)] 
    133194; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite. 
    134195; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite 
    135    points1 = complex(glam[   double(points)], gphi[   double(points)]) 
    136    points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) 
     196  points1 = complex(glam[   double(points)], gphi[   double(points)]) 
     197  points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) 
    137198; droites les equations des droites dont on cherche l''intersection 
    138199; avec la section. 
    139    droites = linearequation(points1, points2) 
    140    inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) 
     200  droites = linearequation(points1, points2) 
     201  inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) 
    141202 
    142203; les ccordonnes geographiques des points que l''on cherche sur la section. 
    143    glamaxe = float(inter) 
    144    gphiaxe = imaginary(inter) 
     204  glamaxe = float(inter) 
     205  gphiaxe = imaginary(inter) 
    145206; on les range ds l''ordre croissant entre les bornes de la section 
    146    if strpos(type, 'x') NE -1 then BEGIN  
    147       sort = sort(glamaxe) 
    148       glamaxe = glamaxe[sort] 
    149       inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) 
    150       glamaxe = glamaxe[inbox] 
    151       sort = sort[inbox] 
    152       gphiaxe = gphiaxe[sort] 
    153    ENDIF ELSE BEGIN 
    154       sort = sort(gphiaxe) 
    155       gphiaxe = gphiaxe[sort] 
    156       inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) 
    157       gphiaxe = gphiaxe[inbox] 
    158       sort = sort[inbox] 
    159       glamaxe = glamaxe[sort] 
    160    ENDELSE 
    161    points = points[sort] 
    162    points1 = points1[sort] 
    163    points2 = points2[sort] 
    164    inter = inter[sort] 
    165    poids = abs(points2-inter)/abs(points2-points1) 
    166  
    167 ;--------------------- 
    168    array = fitintobox(array) 
    169    if array[0] EQ -1 THEN BEGIN 
    170       res = -1 
    171       return 
    172    ENDIF 
    173    if n_elements(valmask) EQ 0 THEN valmask = 1e20 
    174    taille=size(array) 
    175    if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN 
    176       direc = 't' 
    177       array = grossemoyenne(array, 't') 
    178       taille=size(array) 
    179       jpt = 1 
    180    ENDIF 
    181    case 1 of 
     207  if strpos(type, 'x') NE -1 then BEGIN  
     208    sort = sort(glamaxe) 
     209    glamaxe = glamaxe[sort] 
     210    inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) 
     211    glamaxe = glamaxe[inbox] 
     212    sort = sort[inbox] 
     213    gphiaxe = gphiaxe[sort] 
     214  ENDIF ELSE BEGIN 
     215    sort = sort(gphiaxe) 
     216    gphiaxe = gphiaxe[sort] 
     217    inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) 
     218    gphiaxe = gphiaxe[inbox] 
     219    sort = sort[inbox] 
     220    glamaxe = glamaxe[sort] 
     221  ENDELSE 
     222  points = points[sort] 
     223  points1 = points1[sort] 
     224  points2 = points2[sort] 
     225  inter = inter[sort] 
     226  poids = abs(points2-inter)/abs(points2-points1) 
     227 
     228;--------------------- 
     229  array = litchamp(field) 
     230  array = fitintobox(array) 
     231  if array[0] EQ -1 THEN BEGIN 
     232    res = -1 
     233    return 
     234  ENDIF 
     235  if n_elements(valmask) EQ 0 THEN valmask = 1e20 
     236  taille = size(array) 
     237  if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN 
     238    direc = 't' 
     239    array = grossemoyenne(array, 't') 
     240    taille = size(array) 
     241    jpt = 1 
     242  ENDIF 
     243  case 1 of 
    182244;---------------------------------------------------------------------------- 
    183245;xy 
    184246;---------------------------------------------------------------------------- 
    185       taille[0] EQ 2:BEGIN 
    186          value1 = array[double(points)] 
    187          terre = where(value1 GT valmask/10) 
    188          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    189          value2 = array[imaginary(points)] 
    190          terre = where(value2 GT valmask/10) 
    191          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    192          res = poids*value1+(1-poids)*value2 
    193       END 
     247    taille[0] EQ 2:BEGIN 
     248      value1 = array[double(points)] 
     249      terre = where(value1 GT valmask/10) 
     250      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     251      value2 = array[imaginary(points)] 
     252      terre = where(value2 GT valmask/10) 
     253      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     254      res = poids*value1+(1-poids)*value2 
     255    END 
    194256;---------------------------------------------------------------------------- 
    195257;xyz 
    196258;---------------------------------------------------------------------------- 
    197       taille[0] EQ 3 AND jpt EQ 1:BEGIN 
    198          npoints = n_elements(points)  
    199          index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
    200          value1 = array[index] 
    201          terre = where(value1 GT valmask/10) 
    202          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    203          index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
    204          value2 = array[index] 
    205          terre = where(value2 GT valmask/10) 
    206          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    207          poids = poids#replicate(1, nz) 
    208          res = poids*value1+(1-poids)*value2 
     259    taille[0] EQ 3 AND jpt EQ 1:BEGIN 
     260      npoints = n_elements(points)  
     261      index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
     262      value1 = array[index] 
     263      terre = where(value1 GT valmask/10) 
     264      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     265      index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) 
     266      value2 = array[index] 
     267      terre = where(value2 GT valmask/10) 
     268      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     269      poids = poids#replicate(1, nz) 
     270      res = poids*value1+(1-poids)*value2 
    209271; moyenne suivant z ? 
    210          if strpos(type, 'z') EQ -1 then begin 
    211             nan = where(finite(res) EQ 0) 
    212             if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] 
    213             weight = replicate(1, npoints)#e3 
    214             if nan[0] NE -1 then weight[nan] = !values.f_nan 
    215             totalweight = total(weight, 2, /nan) 
    216             zero = where(totalweight EQ 0) 
    217             if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
    218             res = total(res*weight, 2, /nan)/totalweight 
    219             direc = 'z'+string(byte(testvar(var=toto))) 
    220          endif 
    221       END 
     272      if strpos(type, 'z') EQ -1 then begin 
     273        nan = where(finite(res) EQ 0) 
     274        if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt] 
     275        weight = replicate(1, npoints)#e3 
     276        if nan[0] NE -1 then weight[nan] = !values.f_nan 
     277        totalweight = total(weight, 2, /nan) 
     278        zero = where(totalweight EQ 0) 
     279        if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
     280        res = total(res*weight, 2, /nan)/totalweight 
     281        direc = 'z'+string(byte(testvar(var = toto))) 
     282      endif 
     283    END 
    222284;---------------------------------------------------------------------------- 
    223285;xyt 
    224286;---------------------------------------------------------------------------- 
    225       taille[0] EQ 3 AND jpt NE 1:BEGIN 
    226          npoints = n_elements(points)  
    227          index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
    228          value1 = array[index] 
    229          terre = where(value1 GT valmask/10) 
    230          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    231          index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
    232          value2 = array[index] 
    233          terre = where(value2 GT valmask/10) 
    234          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    235          poids = poids#replicate(1, jpt) 
    236          res = poids*value1+(1-poids)*value2 
    237       END 
     287    taille[0] EQ 3 AND jpt NE 1:BEGIN 
     288      npoints = n_elements(points)  
     289      index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
     290      value1 = array[index] 
     291      terre = where(value1 GT valmask/10) 
     292      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     293      index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) 
     294      value2 = array[index] 
     295      terre = where(value2 GT valmask/10) 
     296      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     297      poids = poids#replicate(1, jpt) 
     298      res = poids*value1+(1-poids)*value2 
     299    END 
    238300;---------------------------------------------------------------------------- 
    239301;xyzt 
    240302;---------------------------------------------------------------------------- 
    241       taille[0] EQ 4:BEGIN 
    242          npoints = n_elements(points)  
    243          index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
    244          index = reform(index, npoints, nz, jpt, /over) 
    245          value1 = array[index] 
    246          terre = where(value1 GT valmask/10) 
    247          if terre[0] NE -1 then value1[terre] = !values.f_nan 
    248          index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
    249          index = reform(index, npoints, nz, jpt, /over) 
    250          value2 = array[index] 
    251          terre = where(value2 GT valmask/10) 
    252          if terre[0] NE -1 then value2[terre] = !values.f_nan 
    253          poids = poids#replicate(1, nz*jpt) 
    254          poids = reform(poids, npoints, nz, jpt, /over) 
    255          res = poids*value1+(1-poids)*value2 
     303    taille[0] EQ 4:BEGIN 
     304      npoints = n_elements(points)  
     305      index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
     306      index = reform(index, npoints, nz, jpt, /over) 
     307      value1 = array[index] 
     308      terre = where(value1 GT valmask/10) 
     309      if terre[0] NE -1 then value1[terre] = !values.f_nan 
     310      index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) 
     311      index = reform(index, npoints, nz, jpt, /over) 
     312      value2 = array[index] 
     313      terre = where(value2 GT valmask/10) 
     314      if terre[0] NE -1 then value2[terre] = !values.f_nan 
     315      poids = poids#replicate(1, nz*jpt) 
     316      poids = reform(poids, npoints, nz, jpt, /over) 
     317      res = poids*value1+(1-poids)*value2 
    256318; moyenne suivant z ? 
    257          if strpos(type, 'z') EQ -1 then begin 
    258             nan = where(finite(res) EQ 0) 
    259             if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] 
    260             weight = replicate(1, npoints)#e3 
    261             weight = weight[*]#replicate(1, jpt) 
    262             weight = reform(weight, npoints, nz, jpt, /over) 
    263             if nan[0] NE -1 then weight[nan] = !values.f_nan 
    264             totalweight = total(weight, 2, /nan) 
    265             zero = where(totalweight EQ 0) 
    266             if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
    267             res = total(res*weight, 2, /nan)/totalweight 
    268             direc = 'z'+string(byte(testvar(var=toto))) 
    269          endif 
    270       END 
    271    endcase 
    272 ;--------------------- 
    273  
    274    terre = where(finite(res) EQ 0) 
    275    if terre[0] NE -1 then res[terre] = valmask 
    276  
    277  
    278 ;    plt, tmask[*,*,0],/nocouleur, /rempli, title = '', subtitle = '' 
    279 ;    !p.title='' 
    280 ;    !p.subtitle='' 
    281 ;    dessinetri 
    282 ;    plot,[endpoints[0],endpoints[2]],[endpoints[1],endpoints[3]],/noerase,color=50 
    283  
    284 ;    plot,float(points1),imaginary(points1),/noerase,color=150,psym=1 
    285 ;    plot,float(points2),imaginary(points2),/noerase,color=150,psym=1 
    286 ;    plot,float(inter),imaginary(inter),/noerase,color=250,psym=1 
    287  
    288 ;    plot,[float(points1[0])],[imaginary(points1[0])],/noerase,color=150,psym=1 
    289 ;    plot,[float(points2[0])],[imaginary(points2[0])],/noerase,color=150,psym=1 
    290 ;    plot,[float(inter[0])],[imaginary(inter[0])],/noerase,color=250,psym=1 
    291  
    292  
    293 ;--------------------- 
    294    return 
     319      if strpos(type, 'z') EQ -1 then begin 
     320        nan = where(finite(res) EQ 0) 
     321        if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt] 
     322        weight = replicate(1, npoints)#e3 
     323        weight = weight[*]#replicate(1, jpt) 
     324        weight = reform(weight, npoints, nz, jpt, /over) 
     325        if nan[0] NE -1 then weight[nan] = !values.f_nan 
     326        totalweight = total(weight, 2, /nan) 
     327        zero = where(totalweight EQ 0) 
     328        if zero[0] NE -1 then totalweight[zero] = !values.f_nan 
     329        res = total(res*weight, 2, /nan)/totalweight 
     330        direc = 'z'+string(byte(testvar(var = toto))) 
     331      endif 
     332    END 
     333  endcase 
     334;--------------------- 
     335 
     336  terre = where(finite(res) EQ 0) 
     337  if terre[0] NE -1 then res[terre] = valmask 
     338 
     339  if n_elements(showbuild) then BEGIN  
     340    winsave = !window 
     341    psave = !p 
     342    xsave = !x 
     343    ysave = !y 
     344    plt, findgen(nx, ny), /nodata, /nofill, /rempli, title = '', subtitle = '', coast_thick = 2, window = showbuild 
     345    !p.title = '' 
     346    !p.subtitle = '' 
     347 
     348    plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50 
     349    plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50, psym = 2, thick = 2 
     350 
     351    FOR i = 0, n_elements(points1)-1 DO $ 
     352      plots, [float(points1[i]), float(points2[i])] $ 
     353             , [imaginary(points1[i]), imaginary(points2[i])], color = 150 
     354 
     355    plots, float(points1), imaginary(points1), color = 150, psym = 1 
     356    plots, float(points2), imaginary(points2), color = 150, psym = 1 
     357    plots, float(inter), imaginary(inter), color = 250, psym = 1 
     358     
     359    IF terre[0] NE -1 THEN plots, float(inter[terre]), imaginary(inter[terre]), color = 0, psym = 1 
     360     
     361;      dummy = '' 
     362;      read, dummy,  prompt = 'press return to continue' 
     363    IF !d.name EQ 'PS' THEN erase ELSE wset, winsave 
     364    !p = psave 
     365    !x = xsave 
     366    !y = ysave 
     367  ENDIF  
     368 
     369  restoreboxparam, 'boxparam4section.dat' 
     370 
     371;--------------------- 
     372  return 
    295373end 

  • trunk/ToBeReviewed/TRIANGULATION/tracecote.pro

    r27 r29  
    1515; KEYWORD PARAMETERS:  
    1616; 
    17 ;        CONT_THICK: l''epaisseur du trait pour tracer les 
     17;        COAST_COLOR: the color of the coastline. 
     18;                     defaut value is 0 => black 
     19; 
     20;        COAST_THICK: l''epaisseur du trait pour tracer les 
    1821;        continents. par defaut c''est 1. 
     22; 
     23;        /SURFACE_COASTLINE: to draw the furface coast line instead of 
     24;        the coast line at level firstz[tw]. Usefull only for deep 
     25;        plots! 
    1926; 
    2027;        XSEUIL: pour eliminer les segments de cote qui sont trop 
    … …  
    4552;------------------------------------------------------------ 
    4653;------------------------------------------------------------ 
    47 ; PRO tracecote, CONT_THICK = cont_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 
    48 PRO tracecote, _extra = ex 
    49 @common 
     54PRO tracecote, SURFACE_COASTLINE = surface_coastline, _EXTRA = ex 
     55;-------------------------------------------------------------- 
     56; include commons 
     57@cm_4data 
     58@cm_4mesh 
     59  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     60@updatenew 
     61  ENDIF 
     62;-------------------------------------------------------------- 
    5063   tempsun = systime(1)         ; pour key_performance 
    5164   if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
    … …  
    5669;    
    5770   tempdeux = systime(1)        ; pour key_performance =2 
    58    premierx = 0 > (min([premierxt, premierxu])-1) 
    59    dernierx = (max([dernierxt, dernierxu])+1) < (jpi-1) 
    60    premiery = 0 > (min([premieryt, premieryv])-1) 
    61    derniery = (max([dernieryt, dernieryv])+1) < (jpj-1) 
    62    nx = dernierx-premierx+1 
    63    ny = derniery-premiery+1 
     71   firstx = 0 > (min([firstxt, firstxf])-1) 
     72   lastx = (max([lastxt, lastxf])+1) < (jpi-1) 
     73   firsty = 0 > (min([firstyt, firstyf])-1) 
     74   lasty = (max([lastyt, lastyf])+1) < (jpj-1) 
     75   nx = lastx-firstx+1 
     76   ny = lasty-firsty+1 
    6477; quel niveau vertical choisir ? 
    65    IF strupcase(vargrid) eq 'W' THEN BEGIN  
    66       nz = nzw & premierz = premierzw & dernierz = dernierzw 
    67    ENDIF ELSE BEGIN  
    68       nz = nzt & premierz = premierzt & dernierz = dernierzt 
    69    ENDELSE  
    70    if nz eq jpk then nivz = niveau-1 ELSE nivz = nz-1 
     78   IF keyword_set(surface_coastline) THEN firstz = 0 ELSE $ 
     79     IF strupcase(vargrid) eq 'W' THEN firstz = firstzw ELSE firstz = firstzt 
    7180; attribution du masque et des coordonnes delimitant les limites de la 
    7281; terre (coordonnees f) 
    73    mask = tmask[premierx:dernierx, premiery:derniery, niveau-1] 
    74    xf = glamf[premierx:dernierx, premiery:derniery] 
    75    yf = gphif[premierx:dernierx, premiery:derniery] ; 
     82   mask = tmask[firstx:lastx, firsty:lasty, firstz] 
     83   xf = glamf[firstx:lastx, firsty:lasty] 
     84   yf = gphif[firstx:lastx, firsty:lasty] ; 
    7685   IF testvar(var = key_performance) EQ 2 THEN $ 
    7786    print, 'temps tracecote: determiner mask xf yf', systime(1)-tempdeux 
    … …  
    107116   tempvar = SIZE(TEMPORARY(ind)) ; on efface ind 
    108117; 
    109 ; pour eviter d'avoir des traits de cotes qui travessent tout l''ecran 
    110 ; (qui peuvent etre tres proches sur la sphere mais tres eloignes sur 
    111 ; le dessin suivant la projection) on va selectionner uniquements les 
    112 ; segments de trait de cote qui ne sont pas trop grand, cad qui ne 
    113 ; depassent pas la taille de la fenetre/seuil et pour lesquel le 
    114 ; maskque change de valeur entre les extremites du segment. 
    115 ; 
    116 ; on trace les segments verticaux: 
    117 ; 
    118 ;    if NOT keyword_set(xseuil) then xseuil = 5. 
    119 ;    distanceseuil = (!p.position[2]-!p.position[0])/xseuil 
    120 ; ; liste: liste des points i pourlesquels on va tracer un segment entre 
    121 ; ; le point i,j-1 et i,j 
    122 ;    tempdeux = systime(1)        ; pour key_performance =2 
    123 ;    liste = where(mask+shift(mask, -1, 0) EQ 1 $ 
    124 ;                  AND (xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2 LE distanceseuil) 
    125 ;    IF liste[0] NE -1 THEN BEGIN 
    126 ; ; on recupere lx et ly qui sont les indices ds un tableau 2d des 
    127 ; ; points donnes par liste 
    128 ;       ly = liste/nx & lx = temporary(liste)-nx*ly 
    129 ;       indice = where(ly NE 0)   ; on ne prend pas les points concernant 
    130 ;       if indice[0] NE -1 then begin 
    131 ; ; la premiere ligne car ds ce cas le pt j-1 n''est pas definit 
    132 ;          lx = lx[indice] & ly = ly[temporary(indice)] 
    133 ; ; boucle sur les points concernes et trace du segment 
    134 ; ; rq: on utilise plost au lieu de plot car plots est bcp plus rapide. 
    135 ;          IF testvar(var = key_performance) EQ 2 THEN $ 
    136 ;           print, 'temps tracecote: determiner liste des points concernes par un trait vertical', systime(1)-tempdeux 
    137 ;          tempdeux = systime(1)  ; pour key_performance =2 
    138 ;          for pt = 0, n_elements(lx)-1 do BEGIN  
    139 ;             i = lx[pt] & j = ly[pt] 
    140 ;             plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $ 
    141 ;              , color=c_cote,thick=cont_thick, /normal, _extra = ex 
    142 ;          endfor 
    143 ;          IF testvar(var = key_performance) EQ 2 THEN $ 
    144 ;           print, 'temps tracecote: trace des traits verticaux', systime(1)-tempdeux 
    145 ;       endif 
    146 ;    ENDIF 
    147 ; ; 
    148 ; ; pour le trace des segments horizontaux, c''est la meme chose sauf 
    149 ; ; qu'il faut faire attention si on est periodique:  
    150 ; ; 
    151 ; ; si on est periodique on duplique la premiere colonne et on la met a 
    152 ; ; la fin. (ceci est fait non pas pour le shift qui est par defaut 
    153 ; ; periodique mais pour le plots  
    154 ;    tempdeux = systime(1)        ; pour key_performance =2 
    155 ;    if keyword_set(key_periodique) AND nx EQ jpi then begin 
    156 ;       mask = [mask, mask[0, *]] 
    157 ;       xf = [xf, xf[0, *]] 
    158 ;       yf = [yf, yf[0, *]] 
    159 ;       nx = nx+1 
    160 ;    ENDIF 
    161 ;    if NOT keyword_set(yseuil) then yseuil = 5. 
    162 ;    distanceseuil = (!p.position[3]-!p.position[1])/yseuil 
    163 ;    liste = where(mask+shift(mask, 0, -1) EQ 1 $ 
    164 ;                  AND (xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2 LE distanceseuil) 
    165 ;    IF liste[0] NE -1 THEN BEGIN 
    166 ;       ly = liste/nx & lx = temporary(liste)-nx*ly 
    167 ;       indice = where(ly NE ny-1 AND lx NE 0) 
    168 ;       if indice[0] NE -1 then begin 
    169 ; ; on ne prend pas les points de la  
    170 ; ; premiere colonne et de la derniere ligne (car on l''a rajoute artificiellement!)) 
    171 ;          lx = lx[indice] & ly = ly[temporary(indice)] 
    172 ;          IF testvar(var = key_performance) EQ 2 THEN $ 
    173 ;           print, 'temps tracecote: determiner liste des points concernes par un trait horizontal', systime(1)-tempdeux 
    174 ;          tempdeux = systime(1)  ; pour key_performance =2 
    175 ;          for pt = 0, n_elements(lx)-1 do BEGIN  
    176 ;             i = lx[pt] & j = ly[pt] 
    177 ;             plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $ 
    178 ;              , color=c_cote,thick=cont_thick, /normal, _extra = ex 
    179 ;          endfor 
    180 ;          IF testvar(var = key_performance) EQ 2 THEN $ 
    181 ;           print, 'temps tracecote: trace des traits horizontaux', systime(1)-tempdeux 
    182 ;       endif 
    183 ;    endif 
    184118   if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
    185119   case key_gridtype of 

  • trunk/ToBeReviewed/TRIANGULATION/tracemask.pro

    r27 r29  
    1717; KEYWORD PARAMETERS: 
    1818; 
    19 ;        CONT_THICK: l''epaisseur du trait pour tracer les 
     19;        COAST_COLOR: the color of the coastline. 
     20;                     defaut value is 0 => black 
     21; 
     22;        COAST_THICK: l''epaisseur du trait pour tracer les 
    2023;        continents. par defaut c''est 1. 
    2124; 
    … …  
    3639;------------------------------------------------------------ 
    3740;------------------------------------------------------------ 
    38 PRO tracemask, maskentree, xin, yin, CONT_THICK = cont_thick, OVERPLOT = overplot, _extra = ex 
    39    xentree = xin 
    40    yentree = yin 
     41PRO tracemask, maskentree, xin, yin, COAST_COLOR = coast_color, COAST_THICK = coast_thick, OVERPLOT = overplot, _extra = ex 
     42; 
    4143   if keyword_set(overplot) then return 
    42 @common 
     44;--------------------------------------------------------- 
     45@cm_general 
     46  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     47@updatekwd 
     48  ENDIF 
     49;---------------------------------------------------------  
    4350   tempsun = systime(1)         ; pour key_performance 
    4451; on s''afranchit des problemes de bord: 
    … …  
    4754   nx = tailleentree[1]+1 
    4855   ny = tailleentree[2]+1 
    49 ; on agrandi le mask de une colonne a gauche et de une colonne a droite. 
     56; we check the input axis 
     57  IF n_elements(xin) EQ 0 THEN xentree = findgen(nx-1) ELSE xentree = xin 
     58  IF (size(xentree))[0] EQ 1 THEN xentree = xentree#replicate(1,ny-1) 
     59  IF n_elements(yin) EQ 0 THEN yentree = findgen(ny-1) ELSE yentree = yin 
     60  IF (size(yentree))[0] EQ 1 THEN yentree = replicate(1,nx-1)#yentree 
     61; on agrandi le mask de une colonne a gauche et de une colonne en bas 
    5062   mask = intarr(tailleentree[1]+1, tailleentree[2]+1) 
    5163   mask[1:tailleentree[1], 1:tailleentree[2]] = maskentree 
    … …  
    105117; boucle sur les points concernes et trace du segment 
    106118; rq: on utilise plots au lieu de plot car plots est bcp plus rapide. 
    107          for pt = 0, n_elements(lx)-1 do BEGIN  
     119         for pt = 0L, n_elements(lx)-1 do BEGIN  
    108120            i = lx[pt] & j = ly[pt] 
    109121            plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $ 
    110              , color=c_cote,thick=cont_thick, _extra = ex 
     122              , color = coast_color, thick = coast_thick, _extra = ex 
    111123            if pt LT 5 then begin 
    112124            endif 
    … …  
    129141       print, 'temps tracemask: liste traits horizontaux', systime(1)-tempdeux 
    130142      tempdeux = systime(1)     ; pour key_performance =2 
    131       for pt = 0, n_elements(lx)-1 do BEGIN  
     143      for pt = 0L, n_elements(lx)-1 do BEGIN  
    132144         i = lx[pt] & j = ly[pt] 
    133145         plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $ 
    134           , color=c_cote,thick=cont_thick, _extra = ex 
     146           , color = coast_color, thick = coast_thick, _extra = ex 
    135147      endfor 
    136148      IF testvar(var = key_performance) EQ 2 THEN $ 

  • trunk/ToBeReviewed/TRIANGULATION/triangule.pro

    r27 r29  
    1 FUNCTION triangule, maskentree, REGULIER = regulier, _extra = ex 
     1FUNCTION triangule, maskentree, BASIC = basic, COINMONTE = coinmonte, COINDESCEND = coindescend, _extra = ex 
    22@common 
    3    if keyword_set(regulier) then return, triangule_c(maskentree, REGULIER = regulier, _extra = ex) 
     3; 
     4  IF jpi EQ 1 OR jpj EQ 1 THEN return, -1 
     5; 
     6  IF arg_present(coinmonte) THEN coinmonte = 1 
     7  IF arg_present(coindescend) THEN coindescend = 1 
     8; 
     9   if keyword_set(basic) then $ 
     10     return, triangule_c(maskentree, /BASIC, COINMONTE = coinmonte $ 
     11                         , COINDESCEND = coindescend, _extra = ex) 
     12; 
    413   if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 
    514   if n_elements(maskentree) EQ 0 then maskentree = tmask[*, *, 0] 
    615   case key_gridtype of 
    716      'e':res = triangule_e(maskentree, _extra = ex) 
    8       'c':res = triangule_c(maskentree, REGULIER = regulier, _extra = ex) 
     17      'c':res = triangule_c(maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, _extra = ex) 
    918   endcase 
    1019   return, res 

  • trunk/ToBeReviewed/TRIANGULATION/triangule_c.pro

    r27 r29  
    4242; 
    4343; KEYWORD PARAMETERS: 
    44 ;       /REGULIER: specifie que le masque est sur une grille reguliere 
     44; 
     45;       /BASIC: specifie que le masque est sur une grille basice 
    4546;       (utiliser pour la triangulation ds les coupes verticales et 
    4647;       des hovmoellers) 
     48; 
     49;       /KEEP_CONT: to keep the triangulation even on the continents 
    4750; 
    4851;       COINMONTE=tableau, pour obtenir le tableau de "coins de terre 
    4952;       montant" a traiter avec completecointerre.pro ds la variable 
    5053;       tableau plutot que de la faire passer par la variable globale 
    51 ;       cointerremont. 
     54;       twin_corners_up. 
    5255; 
    5356;       COINDESCEND=tableau cf COINMONTE 
    … …  
    7881;------------------------------------------------------------ 
    7982;------------------------------------------------------------ 
    80 FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $ 
    81                   , REGULIER = regulier, PERIODIQUE = periodique 
     83FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, BASIC = basic, KEEP_CONT = keep_cont 
    8284   tempsun = systime(1)         ; pour key_performance 
    83 @common 
    84    ;; 
     85;--------------------------------------------------------- 
     86@cm_4mesh 
     87  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     88@updatenew 
     89  ENDIF 
    8590;------------------------------------------------------------ 
    8691; le masque est donne ou il faut prendre tmask? 
    … …  
    9297   ny = taille[2] 
    9398; 
    94 ;------------------------------------------------------------ 
    95    if n_elements(periodique) EQ 0 then periodique = keyword_set(key_periodique) 
    96    if keyword_set(key_periodique) and keyword_set(periodique) $ 
    97     AND NOT keyword_set(regulier) then BEGIN  
     99   IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular 
     100;------------------------------------------------------------ 
     101   if keyword_set(key_periodic)*(nx EQ jpi) $ 
     102    AND NOT keyword_set(basic) then BEGIN  
    98103      msk = [msk, msk[0, *]] 
    99104      nx = nx+1 
    … …  
    138143    print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux 
    139144; 
    140    if NOT keyword_set(regulier) then begin 
     145   if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 
    141146      tempdeux = systime(1)     ; pour key_performance =2 
    142147;2 points terre en diagonale montante avec 2 points mer sur la diagonale descendante 
    … …  
    156161       print, 'temps triangule: trouver coindesc', systime(1)-tempdeux 
    157162; 
    158    endif 
     163    ENDIF 
    159164; 
    160165   if n_elements(pts_downward) EQ 1 then BEGIN  
    … …  
    179184      pts_downward =different(pts_downward,derniere_colonne ) 
    180185      pts_downward =different(pts_downward,derniere_ligne ) 
    181       if NOT keyword_set(regulier) then begin 
     186      if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 
    182187         if coinmont[0] NE -1 then begin 
    183188            coinmont =different(coinmont,derniere_colonne ) 
    … …  
    208213;  ca devrait aller mieux dans les prochaines versions d''IDL... 
    209214; 
    210    tempdeux = systime(1)        ; pour key_performance =2 
     215   if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin 
     216      tempdeux = systime(1)     ; pour key_performance =2 
    211217; on enleve les rectangles qui sont entierement dans la terre 
    212    recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 
    213    IF testvar(var = key_performance) EQ 2 THEN $ 
    214     print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 
     218      recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 
     219      IF testvar(var = key_performance) EQ 2 THEN $ 
     220       print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 
    215221 
    216222; en attendant une version qui marche parfaitement, on est contraint 
    … …  
    219225; commun. 
    220226; t1 = systime(1) 
    221    indice = intarr(nx, ny) 
    222    trimask = intarr(nx, ny) 
    223    trimask[0:nx-2, 0:ny-2] = 1 
    224    IF recdsterre[0] NE -1 then BEGIN  
    225       tempdeux = systime(1)     ; pour key_performance =2 
    226       indice[recdsterre] = 1 
    227       if NOT keyword_set(regulier) then begin 
     227      indice = intarr(nx, ny) 
     228      trimask = intarr(nx, ny) 
     229      trimask[0:nx-2, 0:ny-2] = 1 
     230      IF recdsterre[0] NE -1 then BEGIN  
     231         tempdeux = systime(1)  ; pour key_performance =2 
     232         indice[recdsterre] = 1 
     233;      if NOT keyword_set(basic) then begin 
    228234         vire1 = 0 
    229235         vire2 = 0 
    … …  
    247253         IF testvar(var = key_performance) EQ 2 THEN $ 
    248254          print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux 
     255;      endif 
     256         indice[*, ny-1] = 1    ; la deriere colonne te la derniere ligne 
     257         indice[nx-1, *] = 1    ; ne peuvent definir de rectangle. 
     258; 
     259         tempdeux = systime(1)  ; pour key_performance =2 
     260         recgarde = where(indice EQ 0) 
     261; on recupere les numeros des triangles que l'' on va garder 
     262         trigarde = 2*[recgarde-recgarde/nx] 
     263         trigarde = transpose(temporary(trigarde)) 
     264         trigarde = [trigarde, trigarde+1] 
     265;  
     266         triang = triang[*, temporary(trigarde[*])] 
     267         IF testvar(var = key_performance) EQ 2 THEN $ 
     268          print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux 
    249269      endif 
    250       indice[*, ny-1] = 1       ; la deriere colonne te la derniere ligne 
    251       indice[nx-1, *] = 1       ; ne peuvent definirde rectangle. 
    252 ; 
    253       tempdeux = systime(1)     ; pour key_performance =2 
    254       recgarde = where(indice EQ 0) 
    255 ; on recupere les numeros des triangles que l'' on va garder 
    256       trigarde = 2*[recgarde-recgarde/nx] 
    257       trigarde = [trigarde, trigarde+1] 
    258       trigarde = trigarde[sort(trigarde)] 
    259 ;  
    260       triang = triang[*, trigarde] 
    261       IF testvar(var = key_performance) EQ 2 THEN $ 
    262        print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux 
    263270   endif 
    264271; print, 'temps tri triangles', systime(1)-t1  
    265272;------------------------------------------------------------ 
    266 ; quand key_periodique eq 1, triang est une liste d''indice d'un 
     273; quand key_periodic eq 1, triang est une liste d''indice d'un 
    267274; tableau qui a une colonne de trop. 
    268275; il faut ramener ca a la matrice initiale en mettant les indivces de 
    … …  
    270277;------------------------------------------------------------ 
    271278   tempdeux = systime(1)        ; pour key_performance =2 
    272    if keyword_set(key_periodique) and keyword_set(periodique) $ 
    273     AND NOT keyword_set(regulier) then BEGIN  
     279   if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 
     280    AND NOT keyword_set(basic) then BEGIN  
    274281      indicey = triang/nx 
    275282      indicex = triang-indicey*nx 
    … …  
    302309 
    303310;------------------------------------------------------------ 
    304    if arg_present(coinmonte) THEN coinmonte = coinmont ELSE cointerremont = coinmont 
    305    if arg_present(coindescend) THEN coindescend = coindesc ELSE cointerredesc = coindesc 
    306 ; 
    307 ;  
    308 ;------------------------------------------------------------ 
    309  
     311   if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 
     312   if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 
     313;------------------------------------------------------------ 
     314  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     315   @updateold 
     316 ENDIF  
    310317 
    311318   IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun  

  • trunk/ToBeReviewed/TRIANGULATION/triangule_e.pro

    r27 r29  
    3232;------------------------------------------------------------ 
    3333FUNCTION triangule_e, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $ 
    34                   , SHIFTED = shifted, REGULIER = regulier, PERIODIQUE = periodique 
     34                  , SHIFTED = shifted, BASIC = basic 
     35;--------------------------------------------------------- 
     36@cm_4mesh 
     37  IF NOT keyword_set(key_forgetold) THEN BEGIN 
     38@updatenew 
     39  ENDIF 
     40;--------------------------------------------------------- 
    3541   tempsun = systime(1)         ; pour key_performance 
    36 @common 
    37    ;; 
    3842;------------------------------------------------------------ 
    3943; le masque est donne ou il faut prendre tmask? 
    … …  
    4549   ny = sizem[2] 
    4650;------------------------------------------------------------ 
    47    if n_elements(periodique) EQ 0 then periodique = keyword_set(key_periodique) 
    48    if keyword_set(key_periodique) and keyword_set(periodique) $ 
    49     AND NOT keyword_set(regulier) then BEGIN  
     51   if keyword_set(key_periodic)*(nx EQ jpi) $ 
     52    AND NOT keyword_set(basic) then BEGIN  
    5053      msk = [msk, msk[0, *]] 
    5154      nx = nx+1 
    … …  
    132135; ; 
    133136;------------------------------------------------------------ 
    134 ; quand key_periodique eq 1, triang est une liste d''indice d'un 
     137; quand key_periodic eq 1, triang est une liste d''indice d'un 
    135138; tableau qui a une colonne de trop. 
    136139; il faut ramener ca a la matrice initiale en mettant les indivces de 
    … …  
    138141;------------------------------------------------------------ 
    139142   tempdeux = systime(1)        ; pour key_performance =2 
    140    if keyword_set(key_periodique) and keyword_set(periodique) $ 
    141     AND NOT keyword_set(regulier) then BEGIN  
     143   if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 
     144    AND NOT keyword_set(basic) then BEGIN  
    142145      indicey = triang/nx 
    143146      indicex = triang-indicey*nx 
    … …  
    170173 
    171174;------------------------------------------------------------ 
    172 ;    if arg_present(coinmonte) THEN coinmonte = coinmont ELSE cointerremont = coinmont 
    173 ;    if arg_present(coindescend) THEN coindescend = coindesc ELSE cointerredesc = coindesc 
     175;    if arg_present(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 
     176;    if arg_present(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 
    174177; 
    175 ;  
     178;   IF NOT keyword_set(key_forgetold) THEN BEGIN 
     179;    @updateold 
     180;   ENDIF 
    176181;------------------------------------------------------------ 
    177182 
Note: See TracChangeset for help on using the changeset viewer.