source: trunk/ToBeReviewed/PLOTS/VECTEUR/ajoutvect.pro @ 41

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:ajoutvect
;
; PURPOSE:surimprimme des vecteur sur un champ tarce par plt
;
; CATEGORY:grafique
;
; CALLING SEQUENCE:ajoutvect,vecteur 
;
; INPUTS:
;       vecteur: une structure a 2 elements contenant les 2 matrices U
;       et V des valeurs de la composante zonale et meridienne du
;       champ de vecteurs a tracer.
;       par ex:
;       vecteur={matriceu:lec('unsurface'),matricev:lec('vnsurface')}
;       rq:le nom des elements de vecteur n''a aucune importance.
;       vecteur={u:lec('unsurface'),v:lec('vnsurface')} convient aussi
;
; KEYWORD PARAMETERS:
;
;       /UNVECTSUR:un scalaire n on un tableau a 2 elements [n1,n2].
;         dans le premier cas, on tracera un vecteur sur n suivant les
;         x et les y.
;         dans le second cas on tracera un vecteur sur n1 suivant x
;         et un vecteur sur n2 suivant y
;       Rq; pour tracer tous les vecteurs suivant y et 1 sur 2 suivant
;       x mettre unvectsur=[2,1]
;
;       VECTMIN=norme minimun des vecteurs a tracer 
;
;       VECTMAX=norme minimun des vecteurs a tracer 
;
; OUTPUTS:
;
; COMMON BLOCKS:
;       common.pro              
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       10/3/1999
;                       11/6/1999 compatibilite avec NAN et la lecture
;                       des structures.
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
pro ajoutvect,vecteur, vectlegende, UNVECTSUR=unvectsur,VECTMIN=vectmin, VECTMAX=vectmax, _EXTRA = ex
@common
   tempsun = systime(1)         ; pour key_performance
;----------------------------------------------------------------------------
;
   u = litchamp(vecteur.(0))
   u = checkfield(u, 'plt', TYPE = 'xy', /NOQUESTION)
   v = litchamp(vecteur.(1))
   v = checkfield(v, 'plt', TYPE = 'xy', /NOQUESTION)
;-----------------------------------------------------------
;------------------------------------------------------------
; on recupere les eventuelles info sur les champs
;------------------------------------------------------------
   grilleu = litchamp(vecteur.(0), /grid)
   if grilleu EQ '' then grilleu = 'U'
   grillev = litchamp(vecteur.(1), /grid)
   if grillev EQ '' then grillev = 'V'

   IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
   IF grilleu EQ grillev THEN interpolle  = 0 ELSE interpolle = 1
   if keyword_set(inverse) then begin
      rien = u
      u = v
      v = rien
   endif
;------------------------------------------------------------
; on trouve les points que u et v ont en communs
;------------------------------------------------------------
   if interpolle then begin
      indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
      indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
      indicex = inter(indicexu, indicexv)
      indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
      indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
      indicey = inter(indiceyu, indiceyv)
      nx = n_elements(indicex) 
      ny = n_elements(indicey)
      indice2d = lindgen(jpi, jpj)
      indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
      case 1 of
         (size(u))[0] NE 2 OR (size(v))[0] NE 2: return
         (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
          (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
            if nxu NE nx then $
             if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
            IF nxv NE nx THEN $
             if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
            IF nyu NE ny THEN $
             if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
            IF nyv NE ny THEN $
             if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
         END
         (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
          (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
            u = u[indice2d]
            v = v[indice2d]
         END
         ELSE:BEGIN 
            ras = report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
            return
         end
      endcase
;------------------------------------------------------------------
; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee
;------------------------------------------------------------------
      if ny EQ 1 then begin
         u = reform(u, nx, ny)
         v = reform(v, nx, ny)
      endif
;------------------------------------------------------------------
; construction de u et v aux pts T
;-----------------------------------------------------------
      a=u(0,*)
      u=(u+shift(u,1,0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*)=a
      a=v(*,0)
      v=(v+shift(v,0,1))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0)=a
;----------------------------------------------------------------------------
; attribution du mask et des tableau de longitude et latitude
; on recupere la grille complette pour etablir un grand mask etendu ds les 4
; directions pour couvrir les points pour lesquels un pt terre a ete pris en
; compte (faire un petit dessin...)
;----------------------------------------------------------------------------
      vargrid='T'
      msku = (umask())[indice2d+jpi*jpj*firstzt]
      mskv = (vmask())[indice2d+jpi*jpj*firstzt]
      glam = glamt[indice2d]
      gphi = gphit[indice2d]
      if ny EQ 1 then begin
         msku = reform(msku, nx, ny)
         mskv = reform(mskv, nx, ny)
;          glam = reform(glam, nx, ny)
;          gphi = reform(gphi, nx, ny)
      endif
;-----------------------------------------------------------
; on masque u et v le long des cotes (la on l''on ne peut pas calculer
; la moyenne)
;-----------------------------------------------------------
; extention du mask
      u = u*msku*shift(msku,1,0)
      v = v*mskv*shift(mskv,0,1)
   ENDIF ELSE BEGIN 
      u = u*tmask[firstxt:lastxt,firstyt:lastyt,firstzt]
      v = v*tmask[firstxt:lastxt,firstyt:lastyt,firstzt]
      indice2d = lindgen(jpi, jpj)
      indice2d = indice2d[firstxt:lastxt, firstyt:lastyt]
      nx = nxt
      ny = nyt
   endelse
   tabnorme=sqrt(u^2+v^2)
   nan = where(finite(u, /nan) EQ 1)
   if nan[0] NE -1 then u[nan] = 1e5
   nan = where(finite(v, /nan) EQ 1)
   if nan[0] NE -1 then v[nan] = 1e5
   if keyword_set(vectmin) then BEGIN 
      toosmall=where(tabnorme lt vectmin)
      if toosmall[0] NE -1 then begin
         u[toosmall] = 1e5
         v[toosmall] = 1e5
      ENDIF
   endif
   if keyword_set(vectmax) then BEGIN
      toobig=where(tabnorme gt vectmax)
      if toobig[0] NE -1 then begin
         u[toobig] = 1e5
         v[toobig] = 1e5
      ENDIF
   ENDIF
;-----------------------------------------------------------
; remise d''une grande valeur sur tous les points pour lesquelles on ne
; put faire le calcul
;-----------------------------------------------------------
   if interpolle then t2 = msku*shift(msku,1,0)*mskv*shift(mskv,0,1) $
   ELSE t2 = tmask[firstxt:lastxt,firstyt:lastyt,firstzt]
   if NOT keyword_set(key_periodic) OR nx NE jpi then t2(0, *)=0.
   t2(*,0)=0.
   terre=where(t2 eq 0)
   if terre[0] ne -1 then begin
      u(terre)=1e5
      v(terre)=1e5
   ENDIF
;-----------------------------------------------------------
; tracer qu'un vecteur sur
;-----------------------------------------------------------
   if keyword_set(unvectsur) then BEGIN ;
; indx est un vecteur contenant les numero des colonnes a selectionner 
; indy est un vecteur contenant les numero des lignes a selectionner 
      if n_elements(unvectsur) EQ 1 then begin
         indx = where((lindgen(nx) MOD unvectsur[0]) eq 0)
         indy = where((lindgen(ny) MOD unvectsur[0]) eq 0)
      ENDIF ELSE BEGIN 
         indx = where((lindgen(nx) MOD unvectsur[0]) eq 0)
         indy = where((lindgen(ny) MOD unvectsur[1]) eq 0)
      ENDELSE
; a partir de indx et indy on va construire un tableau d''indices 2d
; qui donnera les indices des points intersections des colonnes
; specifiee par indx 
      indicereduit = indx#replicate(1,n_elements(indy))+nx*replicate(1,n_elements(indx))#indy
; on reduit les tableaux qui vont etre passes a vecteur.
      u = u[indicereduit]
      v = v[indicereduit]
      tabnorme = tabnorme[indicereduit]
;
   endif
;-----------------------------------------------------------
; 
;-----------------------------------------------------------
   if keyword_set(inverse) then begin
      rien = u
      u = v
      v = rien
   endif
;-----------------------------------------------------------
; trace des vecteurs
;----------------------------------------------------------
   vecteur, u, v, tabnorme, indice2d, indicereduit, missing=1e5, _extra = ex
;-----------------------------------------------------------
; on complete la legende
;-----------------------------------------------------------
   if terre[0] ne -1 then mini = min(tabnorme[where(t2 eq 1)], max = maxi, /nan) $
   ELSE mini = min(tabnorme, max = maxi, /nan)
   
   if litchamp(vecteur.(0), /u) NE '' then $
    vectlegende = {minmax:[mini, maxi], unite:litchamp(vecteur.(0), /u)} $
   ELSE vectlegende = {minmax:[mini, maxi], unite:varunit}


sortie:
   if keyword_set(key_performance) NE 0 THEN print, 'temps ajoutvect', systime(1)-tempsun 
   return
end