Changeset 432 for trunk/SRC/ToBeReviewed

Timestamp:

05/14/10 12:11:28 (14 years ago)

Author:

smasson

Message:

add VECTNORMCOLOR keyword

File:

• trunk/SRC/ToBeReviewed/PLOTS/VECTEUR/vecteur.pro (modified) (9 diffs)

trunk/SRC/ToBeReviewed/PLOTS/VECTEUR/vecteur.pro

@@ -134 +134 @@
 ; The color of the arrow. Black by default (color 0)
 ;
+; @keyword VECTNORMCOLOR
+; a 2 elements vector used to specify the color of the vector according to
+; its norm. Vector norm will be scaled to color numbers (1 to 254) by using
+; the following rule: VECTNORMCOLOR[0] corresponds to 1 and VECTNORMCOLOR[1] to 254
+;
 ; @keyword VECTTHICK {default=1}
 ; The thick of the arrow.
@@ -169 +174 @@
 ;-
 PRO vecteur, composanteu, composantev, normevecteur, indice2d, reduitindice2d $
-           , CMREF=cmref, MISSING=missing, NORMEREF=normeref $
-           , VECTCOLOR=vectcolor, VECTTHICK=vectthick $
-           , VECTREFPOS=vectrefpos $
-           , VECTREFFORMAT=vectrefformat, NOVECTREF=novectref $
-           , _EXTRA=extra
+             , CMREF = cmref, MISSING = missing, NORMEREF = normeref $
+             , VECTCOLOR = vectcolor, VECTTHICK = vectthick $
+             , VECTREFPOS = vectrefpos, VECTNORMCOLOR = vectnormcolor $
+             , VECTREFFORMAT = vectrefformat, NOVECTREF = novectref $
+             , _EXTRA = extra
 ;
   compile_opt idl2, strictarrsubs
 ;
 @common
-   tempsun = systime(1)         ; For key_performance
-;
-;
-   taille = size(composanteu)
-   nx = taille[1]
-   ny = taille[2]
-   if n_elements(reduitindice2d) EQ 0 then reduitindice2d = lindgen(nx, ny)
-   zu = composanteu
-   zv = composantev
-   norme = normevecteur
-   taille = size(indice2d)
-   nxgd = taille[1]
-   nygd = taille[2]
-;
-   msk = replicate(1, nx, ny)
-   if keyword_set(missing) then terre = where(abs(zu) GE missing/10) ELSE terre = -1
-   if terre[0] NE -1  then BEGIN
-      msk[terre] = 0
-      zu[terre] = 0
-      zv[terre] = 0
-      norme[terre] = 0
-   ENDIF
+  tempsun = systime(1)          ; For key_performance
+;
+;
+  taille = size(composanteu)
+  nx = taille[1]
+  ny = taille[2]
+  if n_elements(reduitindice2d) EQ 0 then reduitindice2d = lindgen(nx, ny)
+  zu = composanteu
+  zv = composantev
+  norme = normevecteur
+  taille = size(indice2d)
+  nxgd = taille[1]
+  nygd = taille[2]
+;
+  msk = replicate(1, nx, ny)
+  if keyword_set(missing) then terre = where(abs(zu) GE missing/10) ELSE terre = -1
+  if terre[0] NE -1  then BEGIN
+    msk[terre] = 0
+    zu[terre] = 0
+    zv[terre] = 0
+    norme[terre] = 0
+  ENDIF
 ;
 ; Stage 1:
@@ -246 +251 @@
 ;
 ; coordinates of the point T (beginning of the arrow) in spherical coordinates.
-   glam = glamt[indice2d[reduitindice2d]]
-   gphi = gphit[indice2d[reduitindice2d]]
+  glam = glamt[indice2d[reduitindice2d]]
+  gphi = gphit[indice2d[reduitindice2d]]
 ;
 ; Coordinates of the point T (beginning of the arrow) in the cartesian reference.
 ; For the sphere, we use a sphere with a radius of 1.
 ;
-   radius = replicate(1,nx*ny)
-   coord_sphe = transpose([ [glam[*]], [gphi[*]], [radius[*]] ])
-   r = cv_coord(from_sphere=coord_sphe,/to_rect,/degrees)
-;
-   x0 = reform(r[0, *], nx, ny)
-   y0 = reform(r[1, *], nx, ny)
-   z0 = reform(r[2, *], nx, ny)
+  radius = replicate(1, nx*ny)
+  coord_sphe = transpose([ [glam[*]], [gphi[*]], [radius[*]] ])
+  r = cv_coord(from_sphere = coord_sphe, /to_rect, /degrees)
+;
+  x0 = reform(r[0, *], nx, ny)
+  y0 = reform(r[1, *], nx, ny)
+  z0 = reform(r[2, *], nx, ny)
 ;
 ; Stage 1, b)
@@ -270 +275 @@
 ;
 ; definition of nu
-   radius = replicate(1,nxgd*nygd)
-   IF finite(glamu[0]*gphiu[0]) NE 0 THEN $
+  radius = replicate(1, nxgd*nygd)
+  IF finite(glamu[0]*gphiu[0]) NE 0 THEN $
      coord_sphe = transpose([ [(glamu[indice2d])[*]], [(gphiu[indice2d])[*]], [radius[*]] ]) $
-   ELSE coord_sphe = transpose([ [(glamf[indice2d])[*]], [(gphit[indice2d])[*]], [radius[*]] ])
-   r = cv_coord(from_sphere=coord_sphe,/to_rect,/degrees)
+  ELSE coord_sphe = transpose([ [(glamf[indice2d])[*]], [(gphit[indice2d])[*]], [radius[*]] ])
+  r = cv_coord(from_sphere = coord_sphe, /to_rect, /degrees)
 ; coordinates of points of the grid u in cartesian.
-   ux = reform(r[0, *], nxgd, nygd)
-   uy = reform(r[1, *], nxgd, nygd)
-   uz = reform(r[2, *], nxgd, nygd)
+  ux = reform(r[0, *], nxgd, nygd)
+  uy = reform(r[1, *], nxgd, nygd)
+  uz = reform(r[2, *], nxgd, nygd)
 ; calculation of nu
-   nux = ux-shift(ux, 1, 0)
-   nuy = uy-shift(uy, 1, 0)
-   nuz = uz-shift(uz, 1, 0)
+  nux = ux-shift(ux, 1, 0)
+  nuy = uy-shift(uy, 1, 0)
+  nuz = uz-shift(uz, 1, 0)
 ; conditions at extremities.
-   if NOT keyword_set(key_periodic) OR nxgd NE jpi then begin
-      nux[0, *] = nux[1, *]
-      nuy[0, *] = nuy[1, *]
-      nuz[0, *] = nuz[1, *]
-   ENDIF
+  if NOT keyword_set(key_periodic) OR nxgd NE jpi then begin
+    nux[0, *] = nux[1, *]
+    nuy[0, *] = nuy[1, *]
+    nuz[0, *] = nuz[1, *]
+  ENDIF
 ; reduction of the grid
-   nux = nux[reduitindice2d]
-   nuy = nuy[reduitindice2d]
-   nuz = nuz[reduitindice2d]
+  nux = nux[reduitindice2d]
+  nuy = nuy[reduitindice2d]
+  nuz = nuz[reduitindice2d]
 ; definition of nv
-   IF finite(glamv[0]*gphiv[0]) NE 0 THEN $
-   coord_sphe = transpose([ [(glamv[indice2d])[*]], [(gphiv[indice2d])[*]], [radius[*]] ]) $
-   ELSE coord_sphe = transpose([ [(glamt[indice2d])[*]], [(gphif[indice2d])[*]], [radius[*]] ])
-   r = cv_coord(from_sphere=coord_sphe,/to_rect,/degrees)
+  IF finite(glamv[0]*gphiv[0]) NE 0 THEN $
+     coord_sphe = transpose([ [(glamv[indice2d])[*]], [(gphiv[indice2d])[*]], [radius[*]] ]) $
+  ELSE coord_sphe = transpose([ [(glamt[indice2d])[*]], [(gphif[indice2d])[*]], [radius[*]] ])
+  r = cv_coord(from_sphere = coord_sphe, /to_rect, /degrees)
 ; coordinates of points of the grid in cartesian.
-   vx = reform(r[0, *], nxgd, nygd)
-   vy = reform(r[1, *], nxgd, nygd)
-   vz = reform(r[2, *], nxgd, nygd)
+  vx = reform(r[0, *], nxgd, nygd)
+  vy = reform(r[1, *], nxgd, nygd)
+  vz = reform(r[2, *], nxgd, nygd)
 ; calcul of nv
-   nvx = vx-shift(vx, 0, 1)
-   nvy = vy-shift(vy, 0, 1)
-   nvz = vz-shift(vz, 0, 1)
+  nvx = vx-shift(vx, 0, 1)
+  nvy = vy-shift(vy, 0, 1)
+  nvz = vz-shift(vz, 0, 1)
 ; conditions at extremities
-   nvx[*, 0] = nvx[*, 1]
-   nvy[*, 0] = nvy[*, 1]
-   nvz[*, 0] = nvz[*, 1]
+  nvx[*, 0] = nvx[*, 1]
+  nvy[*, 0] = nvy[*, 1]
+  nvz[*, 0] = nvz[*, 1]
 ; reduction of the grid
-   nvx = nvx[reduitindice2d]
-   nvy = nvy[reduitindice2d]
-   nvz = nvz[reduitindice2d]
+  nvx = nvx[reduitindice2d]
+  nvy = nvy[reduitindice2d]
+  nvz = nvz[reduitindice2d]
 ;
 ; normalization
 ;
-   normalise, nux, nuy, nuz
-   normalise, nvx, nvy, nvz
+  normalise, nux, nuy, nuz
+  normalise, nvx, nvy, nvz
 ;
 ; Stage 1, c)
@@ -324 +329 @@
 ; coordinates of the vector V in the cartesian reference
 ;
-   direcx = zu*nux + zv*nvx
-   direcy = zu*nuy + zv*nvy
-   direcz = zu*nuz + zv*nvz
+  direcx = zu*nux + zv*nvx
+  direcy = zu*nuy + zv*nvy
+  direcz = zu*nuz + zv*nvz
 ; normalization of the vector V
-   normalise, direcx, direcy, direcz
+  normalise, direcx, direcy, direcz
 ; on divide by 100
-   direcx = direcx/100.
-   direcy = direcy/100.
-   direcz = direcz/100.
+  direcx = direcx/100.
+  direcy = direcy/100.
+  direcz = direcz/100.
 ;
 ; Stage 1, d)
 ; coordinates of the point of the arrow in the cartesian reference.
 
-   x1 = x0 + direcx
-   y1 = y0 + direcy
-   z1 = z0 + direcz
+  x1 = x0 + direcx
+  y1 = y0 + direcy
+  z1 = z0 + direcz
 
 ; coordinates of the point of the arrow in spherical coordinates.
 
-   coord_rect = transpose([ [x1[*]], [y1[*]], [z1[*]] ])
-   r = cv_coord(from_rect=coord_rect,/to_sphere,/degrees)
-   glam1 = reform(r[0, *], nx, ny)
-   gphi1 = reform(r[1, *], nx, ny)
+  coord_rect = transpose([ [x1[*]], [y1[*]], [z1[*]] ])
+  r = cv_coord(from_rect = coord_rect, /to_sphere, /degrees)
+  glam1 = reform(r[0, *], nx, ny)
+  gphi1 = reform(r[1, *], nx, ny)
 
 ;
@@ -355 +360 @@
 ; we modify it
 ;
-   ind = where(glam1 LT !x.range[0] AND glam1+360. LE !x.range[1])
-   if ind[0] NE -1 then glam1[ind] = glam1[ind]+360.
-   ind = where(glam1 GT !x.range[1] AND glam1-360. GE !x.range[0])
-   if ind[0] NE -1 then glam1[ind] = glam1[ind]-360.
-
-   ind = where(glam LT !x.range[0] AND glam+360. LE !x.range[1])
-   if ind[0] NE -1 then glam[ind] = glam[ind]+360.
-   ind = where(glam  GT !x.range[1] AND glam-360. GE !x.range[0])
-   if ind[0] NE -1 then glam[ind] = glam[ind]-360.
+  ind = where(glam1 LT !x.range[0] AND glam1+360. LE !x.range[1])
+  if ind[0] NE -1 then glam1[ind] = glam1[ind]+360.
+  ind = where(glam1 GT !x.range[1] AND glam1-360. GE !x.range[0])
+  if ind[0] NE -1 then glam1[ind] = glam1[ind]-360.
+
+  ind = where(glam LT !x.range[0] AND glam+360. LE !x.range[1])
+  if ind[0] NE -1 then glam[ind] = glam[ind]+360.
+  ind = where(glam  GT !x.range[1] AND glam-360. GE !x.range[0])
+  if ind[0] NE -1 then glam[ind] = glam[ind]-360.
 ;
 ;
 ; Stage 1, e)
 ;
-   r = convert_coord(glam,gphi,/data,/to_normal)
-   x0 = r[0, *]                 ; normal coordinates of the beginning of the array.
-   y0 = r[1, *]                 ;
-
-   r = convert_coord(glam1,gphi1,/data,/to_normal)
-   x1 = r[0, *]                 ; normal coordinates of the ending of the array (Before scaling).
-   y1 = r[1, *]                 ;
+  r = convert_coord(glam, gphi, /data, /to_normal)
+  x0 = r[0, *]                  ; normal coordinates of the beginning of the array.
+  y0 = r[1, *]                  ;
+
+  r = convert_coord(glam1, gphi1, /data, /to_normal)
+  x1 = r[0, *]                  ; normal coordinates of the ending of the array (Before scaling).
+  y1 = r[1, *]                  ;
 ;
 ; tests to avoid that arrows be drawing out of the domain.
 ;
-   out = where(x0 LT !p.position[0] OR x0 GT !p.position[2]  $
-               OR y0 LT !p.position[1] OR y0 GT !p.position[3])
-   if out[0] NE -1 THEN x0[out] = !values.f_nan
+  out = where(x0 LT !p.position[0] OR x0 GT !p.position[2]  $
+              OR y0 LT !p.position[1] OR y0 GT !p.position[3])
+  if out[0] NE -1 THEN x0[out] = !values.f_nan
 ;
 ; Following projections, there may are points at NaN when we pass in normal coordinates.
 ; We delete these points.
 ;
-   nan = finite(x0*y0*x1*y1)
-   number = where(nan EQ 1)
-   x0 = x0[number] & x1 = x1[number]
-   y0 = y0[number] & y1 = y1[number]
-   msk = msk[number]
-   norme = norme[number]
-;
+  nan = finite(x0*y0*x1*y1)
+  number = where(nan EQ 1)
+  x0 = x0[number] & x1 = x1[number]
+  y0 = y0[number] & y1 = y1[number]
+  msk = msk[number]
+  norme = norme[number]
+;
+  points = where(msk EQ 1)
+  IF points[0] NE -1 THEN BEGIN 
+
+  x0 = x0[points] & x1 = x1[points]
+  y0 = y0[points] & y1 = y1[points]
+  msk = msk[points]
+  norme = norme[points]
+
 ; We define the vector direction in the normalize reference.
 ;
-   dirx = x1-x0
-   diry = y1-y0
+    dirx = x1-x0
+    diry = y1-y0
 ;
 ;We pass in polar coordinates to recuperate the angle which wasb the goal of all the first stage!!!
 ;
-
-   dirpol = cv_coord(from_rect = transpose([ [dirx[*]], [diry[*]] ]), /to_polar)
-   dirpol = msk*dirpol[0, *]
+    dirpol = cv_coord(from_rect = transpose([ [dirx[*]], [diry[*]] ]), /to_polar)
+    dirpol = msk*dirpol[0, *]
 ;
 ; Stage 2
@@ -409 +421 @@
 ; Automatic putting at the scale
 ;
-   if NOT keyword_set(cmref) then BEGIN
+    if NOT keyword_set(cmref) then BEGIN
       mipgsz = min(page_size, max = mapgsz)
       sizexfeuille = mipgsz*key_portrait+mapgsz*(1-key_portrait)
@@ -415 +427 @@
       cmref = 5 > floor(sizexfeuille/10.) < 15
       cmref = cmref/10.
-   ENDIF
-   if NOT keyword_set(normeref) then BEGIN
+    ENDIF
+    if NOT keyword_set(normeref) then BEGIN
       value = max(norme)
       puissance10 = 10.^floor(alog10(value))
       normeref = puissance10*floor(value/puissance10)
-   endif
-   cm = 1.*normeref/cmref
+    endif
+    cm = 1.*normeref/cmref
 ;
 ; We modify the array norme to an element having the value cm be represented
-; by a trait of lenght 1 cm on the paper. Norme contain the norme of vectors
+; by a line of lenght 1 cm on the paper. Norme contains the norm of vectors
 ; we want to draw.
 ;
-   norme = 1/(1.*cm)*norme*cv_cm2normal(dirpol)
-;
+; define colors before norme is changed...
+    IF NOT KEYWORD_SET(vectcolor) THEN vectcolor = 0
+    IF keyword_set(vectnormcolor) THEN BEGIN 
+      mp = projsegment([vectnormcolor], [1, 254], /mp)
+      colors = byte(round(mp[0] * norme +  mp[1] ))
+    ENDIF ELSE colors = byte(vectcolor)
+;
+    norme = 1/(1.*cm)*norme*cv_cm2normal(dirpol)
 ;
 ; Stage 3
-; Now that we have the angle and the norme, we recuperate coordinates in
+; Now that we have the angle and the norm, we recuperate coordinates in
 ; rectangular and we draw arrows.
 ;
-   r = cv_coord(from_polar = transpose([ [dirpol[*]], [norme[*]] ]), /to_rect)
-   composantex = r[0, *]
-   composantey = r[1, *]
-;
-   x1 = x0+composantex
-   y1 = y0+composantey
+    r = cv_coord(from_polar = transpose([ [dirpol[*]], [norme[*]] ]), /to_rect)
+    composantex = r[0, *]
+    composantey = r[1, *]
+;
+    x1 = x0+composantex
+    y1 = y0+composantey
 ;
 ; Drawing
 ;
-   if NOT KEYWORD_SET(vectcolor) then vectcolor = 0
-
-   points = where(msk EQ 1)
-   IF points[0] NE -1 THEN arrow, x0[points], y0[points], x1[points], y1[points], /norm $
-    , hsize = -.2, COLOR = vectcolor, THICK = vectthick
+    points = where(msk EQ 1)
+    IF n_elements(colors) GT 1 THEN BEGIN
+      FOR i = 0, n_elements(colors)-1 DO arrow, x0[i], y0[i], x1[i], y1[i], /norm, hsize = -.2, COLOR = colors[i], THICK = vectthick
+    ENDIF ELSE arrow, x0, y0, x1, y1, /norm, hsize = -.2, COLOR = colors, THICK = vectthick
 ;
 ;  Draw an arrow at the right bottom of the drawing as a caption.
 ;
-   if NOT keyword_set(novectref) then BEGIN
+    if NOT keyword_set(novectref) then BEGIN
       dx = cmref*cv_cm2normal(0) ; Length of the vector of reference in normalized coordinates.
       if keyword_set(vectrefformat) then $
-       normelegende = strtrim(string(normeref, format = vectrefformat), 1)+' ' $
+         normelegende = strtrim(string(normeref, format = vectrefformat), 1)+' ' $
       ELSE normelegende = strtrim(normeref, 1)+' '
 ;
       if keyword_set(vectrefpos) then begin
-         r = convert_coord(vectrefpos,/data, /to_normal)
-         x0 = r[0]
-         y0 = r[1]
+        r = convert_coord(vectrefpos, /data, /to_normal)
+        x0 = r[0]
+        y0 = r[1]
       ENDIF ELSE BEGIN
-         x0 = !x.window[1]-dx
-         r = convert_coord(!d.x_ch_size, !d.y_ch_size, /device, /to_normal)
-         dy = 3*r[1]*!p.charsize
-         y0 = !y.window[0]-dy
+        x0 = !x.window[1]-dx
+        r = convert_coord(!d.x_ch_size, !d.y_ch_size, /device, /to_normal)
+        dy = 3*r[1]*!p.charsize
+        y0 = !y.window[0]-dy
       ENDELSE
 
@@ -471 +488 @@
       xyouts, x0, y0, normelegende, /norm, align = 1, charsize = !p.charsize, color = 0
 
-   endif
-;
-;
-
-   if keyword_set(key_performance) NE 0 THEN print, 'temps vecteur', systime(1)-tempsun
+    endif
+  endif
+;
+;
+
+  if keyword_set(key_performance) NE 0 THEN print, 'temps vecteur', systime(1)-tempsun
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-   return
+  return
 END