source: trunk/SRC/ToBeReviewed/TRIANGULATION/drawcoast_c.pro

;+
; @file_comments
;
; @categories
;
; @param MASK
;
; @param XF
;
; @param YF
;
; @param NX
;
; @param NY
;
; @keyword COAST_COLOR {default=0}
; The color of the coastline.
; Default is black (0).
;
; @keyword COAST_THICK {default=1}
; The thick of the trait to trace continents
;
; @keyword XSEUIL {default=5}
; To eliminate segments of coasts which are to big (which link points which can
; be close on the sphere but distant on the drawing). We delete all segments
; whose the size surpass the size of the window following X/XSEUIL.
; But it can be to big if we do a big zoom or a little one for some
; projections... We specify it the keyword thanks to this keyword!
;
; @keyword YSEUIL
; See XSEUIL
;
; @keyword _EXTRA
; Used to pass keywords
;
; @returns
;
; @uses
;
; @restrictions
;
; @examples
;
; @history
;
; @version
;
; @todo
; Seb: remplir le header
;
;-
PRO drawcoast_c, mask, xf, yf, nx, ny $
               , COAST_COLOR=coast_color, COAST_THICK=coast_thick $
               , YSEUIL=yseuil, XSEUIL=xseuil, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;---------------------------------------------------------
   tempsun = systime(1)         ; For key_performance
;---------------------------------------------------------
;
; We trace vertical segments:
;
   if NOT keyword_set(yseuil) then yseuil = 5. < (min([nx, ny])-2)
   distanceseuil = (!p.position[3]-!p.position[1])/yseuil
; list: list of points i for which we will trace a segment between the point i,j-1 and i,j
   tempdeux = systime(1)        ; For key_performance =2
   liste = where((mask+shift(mask, -1, 0)) EQ 1 $
                 AND ((xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2) LE distanceseuil^2)
   IF liste[0] NE -1 THEN BEGIN
; We recuperate lx an dly which are indexes in a 2d array of points given by list
      ly = liste/nx & lx = temporary(liste)-nx*ly
      indice = where(ly NE 0)   ; We do not take points concerning
      if indice[0] NE -1 then begin
; the first line because in this case, the point j-1 is undefined.
         lx = lx[indice] & ly = ly[temporary(indice)]
; Loop on concerned points and drawing of the segment.
; Comment: we use plots instead of plot because plots goes faster.
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps tracecote: determiner liste des points concernes par un trait vertical', systime(1)-tempdeux
         tempdeux = systime(1)  ; For key_performance =2
         for pt = 0L, n_elements(lx)-1 do BEGIN
            i = lx[pt] & j = ly[pt]
            plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $
              , color = coast_color, thick = coast_thick, /normal, _extra = ex
         endfor
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps tracecote: trace des traits verticaux', systime(1)-tempdeux
      endif
   ENDIF
;
; For the drawing of horizontal segments , it is the same thing but we have to be careful if it is periodic.
;
; If it is periodic, we duplicate the first column and we put it at the end.
; (This is made not for the shift which is periodic by default but for the plots)
   tempdeux = systime(1)        ; For key_performance =2
   if keyword_set(key_periodic) AND nx EQ jpi then begin
      mask = [mask, mask[0, *]]
      xf = [xf, xf[0, *]]
      yf = [yf, yf[0, *]]
      nx = nx+1
   ENDIF
   if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2)
   distanceseuil = (!p.position[2]-!p.position[0])/xseuil
   liste = where((mask+shift(mask, 0, -1)) EQ 1 $
                 AND ((xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2) LE distanceseuil^2)
   IF liste[0] NE -1 THEN BEGIN
      ly = liste/nx & lx = temporary(liste)-nx*ly
      indice = where(ly NE ny-1 AND lx NE 0)
      if indice[0] NE -1 then begin
; We do not take points of the first column and of the last line (because we added artificially)
         lx = lx[indice] & ly = ly[temporary(indice)]
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps tracecote: determiner liste des points concernes par un trait horizontal', systime(1)-tempdeux
         tempdeux = systime(1)  ; For key_performance =2
         for pt = 0L, n_elements(lx)-1 do BEGIN
            i = lx[pt] & j = ly[pt]
            plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $
              , color = coast_color, thick = coast_thick, /normal, _extra = ex
         endfor
         IF testvar(var = key_performance) EQ 2 THEN $
          print, 'temps tracecote: trace des traits horizontaux', systime(1)-tempdeux
      endif
   endif
;---------------------------------------------------------
   if keyword_set(key_performance) THEN print, 'temps drawcoast_c', systime(1)-tempsun
   return
end