source: trunk/PLOTS/VECTEUR/vecteur.pro @ 2

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:vecteur
;
; PURPOSE: trace des vecteurs (meme situees sur une grille tordue) sur
; n''importe quelle projection, de telle sorte que tous les vecteurs
; aient une norme comparable sur le dessin (en clair, un vecteur qui
; doit faire 1cm le fait quelque soit la projection et sa position sur
; la sphere) 
;
; CATEGORY:trace de vecteur
;
; CALLING SEQUENCE:vecteur, composanteu, composantev, normevecteur, indice2d, reduitindice2d
; 
; INPUTS: COMPOSANTEU et COMPOSANTEV: ce sont les composantes des
; vecteurs a tracer. Ces tableaux 2d ont la meme dimension que
; reduitindice2d (cf apres).
;       INDICE2D: indice permettant de passer d''un tableau jpi,jpj au
;       zoom surlequel on fait le dessin.
;       REDUITINDICE2D: indice permettant de passer d''un tableau
;       definit par indice2d au tableau pourlequel on a reelement des
;       vecteurs a tracer (en clair: c''est par ex qd on ne trace par exemple
;       que un vecteur sur 2)
;
; KEYWORD PARAMETERS:
;
;      CMREF: la longeur en cm sur le papier que diot faire la fleche
;      de norme normeref. par defaut ajuste au dessin et compris entre
;      .5 et 1.5 cm
;
;      MISSING: la valeur d''une missing value. ne pas utilisder ce
;      mot cle. fixe a 1e5 par ajoutvect.pro
;      
;      NORMEREF: la norme de la fleche de reference. par defaut on
;      essaie de faire qqch qui colle pas trop mal!
;
;      VECTCOLOR: la couleur de la fleche. Par defaut noir (couleur 0)
;
;      VECTTHICK; l''epaissuer de la fleche. par defaut 1.
;
;      VECTREFPOS: vecteur de 2 elements specifiant la position en
;      coordonnees DATA du debut du vecteur de reference. Par defaut
;      en bas a droite du dessin.
;
;      VECTREFFORMAT: format a utiliser pour specifier la norme du
;      vecteur de reference.
;
;      /NOVECTREF: pour supprimer l'affichage du vecteur de reference.
;
; OUTPUTS:
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;  Creation : 13/02/98 G. Roullet (grlod@lodyc.jussieu.fr)
;  Modification : 14/01/99 realise la transformation
;  spherique<->cartesien G. Roullet
;                 12/03/99 verification de la routine G. Roullet
;  8/11/1999:
;  G. Roullet et Sebastien Masson (smasson@lodyc.jussieu.fr)
;  adaptation pour les zoom. reverification...traitement separe de la
;  direction et de la norme des vecteurs. mots cles NORMEREF et CMREF.
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------

FUNCTION cv_cm2normal, angle
;
; donne la longeur en coordonnes normales d''un trait ortiente de
; angle par rapport a l''axe des x et qui doit faire 1 cm sur le
; dessin.
; angle peut etre un tableau.
;
@common
; quelle est la longeur en coordonnees normales d''un trait qui fera 1
; cm sur le parier et qui est parallele a x?
   sizexfeuille = petitfeuille*key_portrait+grandfeuille*(1-key_portrait)
   sizeyfeuille = grandfeuille*key_portrait+petitfeuille*(1-key_portrait)
   cm_en_normal = 1./sizexfeuille
;
;  si le rapport d''aspect de la fenetre n''est pas egale a 1, la
;  longeur en coordonees normalise d''un trait de 1cm varie suivant
;  l''angle polaire de ce trait.
;
   aspect = sizexfeuille/sizeyfeuille
   cm_en_normal = cm_en_normal*sqrt( 1 +(aspect^2-1)*sin(angle)^2 )
;
   return, cm_en_normal
END
;
PRO normalise, u, v, w
;
; normalise le vecteur
;
   IF n_elements(w) NE 0 THEN BEGIN 
      norme = sqrt(u^2.+v^2.+w^2.)
      ind = where(norme NE 0)
      u(ind) = u(ind)/norme[ind]
      v(ind) = v(ind)/norme[ind]
      w(ind) = w(ind)/norme[ind]
   ENDIF ELSE BEGIN
      norme = sqrt(u^2.+v^2.)
      ind = where(norme NE 0)
      u(ind) = u(ind)/norme[ind]
      v(ind) = v(ind)/norme[ind]
   ENDELSE 
END
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
@common
   tempsun = systime(1)         ; pour 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
;
; Etape 1: 
;
; etant donne la direction et le sens que le vecteur a sur la
; sphere, il faut se debrouiller pour determiner cette direction et le
; sens qu'aura le vecteur sur l''ecran ou la feuille une fois qu''il
; aura ete projete.
;
; En theorie: sur la sphere, un vecteur en un point donne a pour
; direction la tangente au cercle qui passe par le centre de la terre
; et par le vecteur. Donc trouver la direction une fois la projection
; faite ds le plan 2d, c''est trouver la tangente a la courbe
; representant la projection du cercle sur le plan 2d au point
; representant la projection du point de depart de la shere sur le
; plan 2d!!!
;
; En pratique on ne connait pas la definition de la courbe que donne
; la projection d''un cercle alors trouver sa tangente en un point...
;
; Ce que l''on fait:
; Ds un repere cartesien 3d, 
;       a) on trouve les coordonnees du point T (debut de la fleche)
;       situe sur la shere.
;       b) pour chaque point T on determine les directions locales
;       definies par la grille en ce point et auxquelles se rapportent
;       les coordonnes (u,v) du vecteur. Ces directions locales sont
;       definies par les gradiants des glam et gphi. Une fois ces
;       directions obtenues, on les considere comme orthogonales et en
;       les normant on construit un repere orthonormal (nu,nv) auquel se
;       rapporte les coordonnes (u,v) du vecteur. Ds le repere
;       cartesien 3d de depart, le vecteur est definit par:
;       V=u*nu+v*nv (ou V, nu et nv sont des vecteurs 3d et u et v des
;       scalaires)
;       c) pour approximer la tangente au cercle par la corde definie
;       par le debut et la fin de la fleche, on va normaliser V puis
;       le diviser par 100.
;       d) ceci nous permet de determiner les coordonnees ds le repere
;       cartesien 3d des extremites de la corde, on les passe en
;       coordonnes sphereriques de facon a recuperer les positions en
;       latitude longitude de ces points sur la sphere.
;       e) on passe les coordonnees de ces points en coordonnees
;       normalise puis en corrdonnes polaire de facon a trouver
;       l''angle et la direction qu''ils determinent sur le dessin...
;
;
; etape 1, a)
;
;
; coordonnes du point T (debut de la fleche) en coordonnes sheriques
   glam = glamt[indice2d[reduitindice2d]]
   gphi = gphit[indice2d[reduitindice2d]]
;
; coordonnes du point T (debut de la fleche) dans le repere cartesien
; on utilise comme shere, une shere de rayon 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)
;
; etape 1, b)
;
; Construction du vecteur nu (resp. nv), vecteur norme porte par
; l''axe des points u(i,j) et u(i-1,j) (resp v(i,j) et v(i,j-1))
; qui definissent pour chaque point sur la shere les directions locales
; associee a u et v. ces vecteurs definissent un repere orthonorme
; local. 
; ces vecteurs sont construits dans un repere cartesien (cv_coord)
; on a choisit un rayon de la terre unite (unit).
;
; definition de nu
   radius = replicate(1,nxgd*nygd)
   coord_sphe = transpose([ [(glamu[indice2d])[*]], [(gphiu[indice2d])[*]], [radius(*)] ])
   r = cv_coord(from_sphere=coord_sphe,/to_rect,/degrees)
; coordonnes de points de la grille u en cartesien
   ux = reform(r(0, *), nxgd, nygd)
   uy = reform(r(1, *), nxgd, nygd)
   uz = reform(r(2, *), nxgd, nygd)
; calcul de nu 
   nux = ux-shift(ux, 1, 0)
   nuy = uy-shift(uy, 1, 0)
   nuz = uz-shift(uz, 1, 0)
; conditions aux limites
   if NOT keyword_set(key_periodique) OR nxgd NE jpi then begin
      nux[0, *] = nux[1, *]
      nuy[0, *] = nuy[1, *]
      nuz[0, *] = nuz[1, *]
   ENDIF
; reduction de la grille
   nux = nux[reduitindice2d]
   nuy = nuy[reduitindice2d]
   nuz = nuz[reduitindice2d]
; definition de nv
   coord_sphe = transpose([ [(glamv[indice2d])[*]], [(gphiv[indice2d])[*]], [radius(*)] ])
   r = cv_coord(from_sphere=coord_sphe,/to_rect,/degrees)
; coordonnes de points de la grille v en cartesien
   vx = reform(r(0, *), nxgd, nygd)
   vy = reform(r(1, *), nxgd, nygd)
   vz = reform(r(2, *), nxgd, nygd)
; calcul de nv 
   nvx = vx-shift(vx, 0, 1)
   nvy = vy-shift(vy, 0, 1)
   nvz = vz-shift(vz, 0, 1)
; conditions aux limites
   nvx[*, 0] = nvx[*, 1]
   nvy[*, 0] = nvy[*, 1]
   nvz[*, 0] = nvz[*, 1]
; reduction de la grille
   nvx = nvx[reduitindice2d]
   nvy = nvy[reduitindice2d]
   nvz = nvz[reduitindice2d]
;
; normalisation
;
   normalise, nux, nuy, nuz
   normalise, nvx, nvy, nvz
;
; etape 1, c)
;
; coordonnes du vecteur V ds le repere cartesion
;
   direcx = zu*nux + zv*nvx
   direcy = zu*nuy + zv*nvy
   direcz = zu*nuz + zv*nvz
; normalisation du vecteur v
   normalise, direcx, direcy, direcz
; on divise par 100.
   direcx = direcx/100.
   direcy = direcy/100.
   direcz = direcz/100.
;
; etape 1, d)
; coordonnees de la pointe de la fleche dans le repere cartesien

   x1 = x0 + direcx
   y1 = y0 + direcy
   z1 = z0 + direcz

; coordonnees de la pointe en spherique

   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)

;
; modif des glam tout se passe bien au niveau de la ligne de
; changement de date...attention, il ne faut pas couper les fleches
; qui sortent de la fenetre!
; test: si il sort du cadre mais qu''avec un +/- 360 il y rentre on le
; modifie...
;
   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.
;
;
; etape 1, e)
;
   r = convert_coord(glam,gphi,/data,/to_normal) 
   x0 = r(0, *)                 ; coordonnes normales du debut de la fleche
   y0 = r(1, *)                 ; 
   
   r = convert_coord(glam1,gphi1,/data,/to_normal) 
   x1 = r(0, *)                 ; coordonnes normales de la fin de la fleche (avant scaling)
   y1 = r(1, *)                 ; 
;
; tests pour eviter que des fleches soient dessineees hors du domaine
;
   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
;
; suivant les projections, il peu y a voir des points a nan qd on
; passe en coordonnes normales. on supprime ces 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]
;
; on definit le vecteur direction dans le repere normalise
;
   dirx = x1-x0
   diry = y1-y0
;
; on passe en polaire pour recuperer l''angle qui etait le but de
; toute la partie 1!!!
;

   dirpol = cv_coord(from_rect = transpose([ [dirx[*]], [diry[*]] ]), /to_polar)
   dirpol = msk*dirpol[0, *]
;
; 2eme etape....
;
; maintenant on s''occupe de la norme...
;
; Mise a l''echelle automatique
;
   if NOT keyword_set(cmref) then BEGIN 
      sizexfeuille = petitfeuille*key_portrait+grandfeuille*(1-key_portrait)
      sizexfeuille = 10.*sizexfeuille
      cmref = 5 > floor(sizexfeuille/10.) < 15
      cmref = cmref/10.
   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
;
; on modifie le tableau norme de facon a ce que un element qui a la
; valeur cm soit represente par un trait de longueur 1cm sur le
; papier.
; norme contient la norme des vecteur que l''on veut dessiner
;
   norme = 1/(1.*cm)*norme*cv_cm2normal(dirpol)
;
;
; 3eme etape....
; maintenant qu''on a l''angle et la norme, et bien on recupere les
; coordonnes en rectangulaire et on dessine les fleches.
;
   r = cv_coord(from_polar = transpose([ [dirpol[*]], [norme[*]] ]), /to_rect)
   composantex = r(0, *)
   composantey = r(1, *)
;
   x1 = x0+composantex
   y1 = y0+composantey
;
; c''est parti pour le trace !
;
   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
;
;  Dessine une fleche en bas a droite de la figure en guise de legende
;
   if NOT keyword_set(novectref) then BEGIN
      dx = cmref*cv_cm2normal(0) ; longuer du vecteur de reference en coordonnes normalisees.
      if keyword_set(vectrefformat) then $
       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]
      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
      ENDELSE

      arrow, x0, y0, x0+dx, y0, /norm, hsize = -.2, color = 0
      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 
;------------------------------------------------------------
;------------------------------------------------------------
   return
END