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

Last change on this file since 2 was 2, checked in by opalod, 22 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 10.0 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:ajoutvect
6;
7; PURPOSE:surimprimme des vecteur sur un champ tarce par plt
8;
9; CATEGORY:grafique
10;
11; CALLING SEQUENCE:ajoutvect,vecteur
12;
13; INPUTS:
14;       vecteur: une structure a 2 elements contenant les 2 matrices U
15;       et V des valeurs de la composante zonale et meridienne du
16;       champ de vecteurs a tracer.
17;       par ex:
18;       vecteur={matriceu:lec('unsurface'),matricev:lec('vnsurface')}
19;       rq:le nom des elements de vecteur n''a aucune importance.
20;       vecteur={u:lec('unsurface'),v:lec('vnsurface')} convient aussi
21;
22; KEYWORD PARAMETERS:
23;
24;       /UNVECTSUR:un scalaire n on un tableau a 2 elements [n1,n2].
25;         dans le premier cas, on tracera un vecteur sur n suivant les
26;         x et les y.
27;         dans le second cas on tracera un vecteur sur n1 suivant x
28;         et un vecteur sur n2 suivant y
29;       Rq; pour tracer tous les vecteurs suivant y et 1 sur 2 suivant
30;       x mettre unvectsur=[2,1]
31;
32;       VECTMIN=norme minimun des vecteurs a tracer
33;
34;       VECTMAX=norme minimun des vecteurs a tracer
35;
36; OUTPUTS:
37;
38; COMMON BLOCKS:
39;       common.pro             
40;
41; SIDE EFFECTS:
42;
43; RESTRICTIONS:
44;
45; EXAMPLE:
46;
47; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
48;                       10/3/1999
49;                       11/6/1999 compatibilite avec NAN et la lecture
50;                       des structures.
51;-
52;------------------------------------------------------------
53;------------------------------------------------------------
54;------------------------------------------------------------
55pro ajoutvect,vecteur, vectlegende, UNVECTSUR=unvectsur,VECTMIN=vectmin, VECTMAX=vectmax, _EXTRA = ex
56@common
57   tempsun = systime(1)         ; pour key_performance
58;----------------------------------------------------------------------------
59;
60   u = litchamp(vecteur.(0))
61   u = checkfield(u, 'plt', TYPE = 'xy', /NOQUESTION)
62   v = litchamp(vecteur.(1))
63   v = checkfield(v, 'plt', TYPE = 'xy', /NOQUESTION)
64;         niv=''
65   IF (size(u))[0] EQ 3 THEN BEGIN
66      niveau = xquestion('Le tableau d''entree est un tableau 3d, a quel niveau faut-il faire les vecteurs? (par defaut le premier niveau) ', strtrim(niveau, 1), /chkwidget)
67      niveau = 1 > fix(niveau) < jpk
68      u=u(*,*,niveau-1)
69      v=v(*,*,niveau-1)
70   endif
71
72;-----------------------------------------------------------
73;------------------------------------------------------------
74; on recupere les eventuelles info sur les champs
75;------------------------------------------------------------
76   grilleu = litchamp(vecteur.(0), /grid)
77   if grilleu EQ '' then grilleu = 'U'
78   grillev = litchamp(vecteur.(1), /grid)
79   if grillev EQ '' then grillev = 'V'
80   IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
81   IF grilleu EQ grillev THEN interpolle  = 0 ELSE interpolle = 1
82   if keyword_set(inverse) then begin
83      rien = u
84      u = v
85      v = rien
86   endif
87;------------------------------------------------------------
88; on trouve les points que u et v ont en communs
89;------------------------------------------------------------
90   if interpolle then begin
91      indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1]
92      indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1]
93      indicex = inter(indicexu, indicexv)
94      indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1]
95      indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1]
96      indicey = inter(indiceyu, indiceyv)
97      nx = n_elements(indicex)
98      ny = n_elements(indicey)
99      indice2d = lindgen(jpi, jpj)
100      indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
101;------------------------------------------------------------
102; extraction de u et v sur le domaine qui convient
103;------------------------------------------------------------
104      case 1 of
105         (size(u))[0] NE 2 OR (size(v))[0] NE 2: return
106         (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
107          (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
108            if nxu NE nx then $
109             if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
110            IF nxv NE nx THEN $
111             if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
112            IF nyu NE ny THEN $
113             if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
114            IF nyv NE ny THEN $
115             if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
116         END
117         (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
118          (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
119            u = u[indice2d]
120            v = v[indice2d]
121         END
122         ELSE:BEGIN
123            ras = report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
124            return
125         end
126      endcase
127;------------------------------------------------------------------
128; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee
129;------------------------------------------------------------------
130      if ny EQ 1 then begin
131         u = reform(u, nx, ny)
132         v = reform(v, nx, ny)
133      endif
134;------------------------------------------------------------------
135; construction de u et v aux pts T
136;-----------------------------------------------------------
137      a=u(0,*)
138      u=(u+shift(u,1,0))/2.
139      if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,*)=a
140      a=v(*,0)
141      v=(v+shift(v,0,1))/2.
142      if NOT keyword_set(key_periodique) OR nx NE jpi then v(*,0)=a
143;----------------------------------------------------------------------------
144; attribution du mask et des tableau de longitude et latitude
145; on recupere la grille complette pour etablir un grand mask etendu ds les 4
146; directions pour couvrir les points pour lesquels un pt terre a ete pris en
147; compte (faire un petit dessin...)
148;----------------------------------------------------------------------------
149      vargrid='T'
150      msku = (umask())[indice2d+jpi*jpj*(niveau-1)]
151      mskv = (vmask())[indice2d+jpi*jpj*(niveau-1)]
152      glam = glamt[indice2d]
153      gphi = gphit[indice2d]
154      if ny EQ 1 then begin
155         msku = reform(msku, nx, ny)
156         mskv = reform(mskv, nx, ny)
157;          glam = reform(glam, nx, ny)
158;          gphi = reform(gphi, nx, ny)
159      endif
160;-----------------------------------------------------------
161; on masque u et v le long des cotes (la on l''on ne peut pas calculer
162; la moyenne)
163;-----------------------------------------------------------
164; extention du mask
165      u = u*msku*shift(msku,1,0)
166      v = v*mskv*shift(mskv,0,1)
167   ENDIF ELSE BEGIN
168      u = u*tmask[premierxt:dernierxt,premieryt:dernieryt,premierzt]
169      v = v*tmask[premierxt:dernierxt,premieryt:dernieryt,premierzt]
170      indice2d = lindgen(jpi, jpj)
171      indice2d = indice2d[premierxt:dernierxt, premieryt:dernieryt]
172      nx = nxt
173      ny = nyt
174   endelse
175   tabnorme=sqrt(u^2+v^2)
176   nan = where(finite(u, /nan) EQ 1)
177   if nan[0] NE -1 then u[nan] = 1e5
178   nan = where(finite(v, /nan) EQ 1)
179   if nan[0] NE -1 then v[nan] = 1e5
180   if keyword_set(vectmin) then BEGIN
181      troppetit=where(tabnorme lt vectmin)
182      if troppetit[0] NE -1 then begin
183         u[troppetit] = 1e5
184         v[troppetit] = 1e5
185      ENDIF
186   endif
187   if keyword_set(vectmax) then BEGIN
188      tropgrand=where(tabnorme gt vectmax)
189      if tropgrand[0] NE -1 then begin
190         u[tropgrand] = 1e5
191         v[tropgrand] = 1e5
192      ENDIF
193   ENDIF
194;-----------------------------------------------------------
195; remise d''une grande valeur sur tous les points pour lesquelles on ne
196; put faire le calcul
197;-----------------------------------------------------------
198   if interpolle then t2 = msku*shift(msku,1,0)*mskv*shift(mskv,0,1) $
199   ELSE t2 = tmask[premierxt:dernierxt,premieryt:dernieryt,premierzt]
200   if NOT keyword_set(key_periodique) OR nx NE jpi then t2(0, *)=0.
201   t2(*,0)=0.
202   terre=where(t2 eq 0)
203   if terre[0] ne -1 then begin
204      u(terre)=1e5
205      v(terre)=1e5
206   ENDIF
207;-----------------------------------------------------------
208; tracer qu'un vecteur sur
209;-----------------------------------------------------------
210   if keyword_set(unvectsur) then BEGIN ;
211; indx est un vecteur contenant les numero des colonnes a selectionner
212; indy est un vecteur contenant les numero des lignes a selectionner
213      if n_elements(unvectsur) EQ 1 then begin
214         indx = where((lindgen(nx) MOD unvectsur[0]) eq 0)
215         indy = where((lindgen(ny) MOD unvectsur[0]) eq 0)
216      ENDIF ELSE BEGIN
217         indx = where((lindgen(nx) MOD unvectsur[0]) eq 0)
218         indy = where((lindgen(ny) MOD unvectsur[1]) eq 0)
219      ENDELSE
220; a partir de indx et indy on va construire un tableau d''indices 2d
221; qui donnera les indices des points intersections des colonnes
222; specifiee par indx
223      indicereduit = indx#replicate(1,n_elements(indy))+nx*replicate(1,n_elements(indx))#indy
224; on reduit les tableaux qui vont etre passes a vecteur.
225      u = u[indicereduit]
226      v = v[indicereduit]
227      tabnorme = tabnorme[indicereduit]
228;
229   endif
230;-----------------------------------------------------------
231;
232;-----------------------------------------------------------
233   if keyword_set(inverse) then begin
234      rien = u
235      u = v
236      v = rien
237   endif
238;-----------------------------------------------------------
239; trace des vecteurs
240;----------------------------------------------------------
241   vecteur, u, v, tabnorme, indice2d, indicereduit, missing=1e5, _extra = ex
242;-----------------------------------------------------------
243; on complete la legende
244;-----------------------------------------------------------
245   if terre[0] ne -1 then mini = min(tabnorme[where(t2 eq 1)], max = maxi, /nan) $
246   ELSE mini = min(tabnorme, max = maxi, /nan)
247   
248   if litchamp(vecteur.(0), /u) NE '' then $
249    vectlegende = {minmax:[mini, maxi], unite:litchamp(vecteur.(0), /u)} $
250   ELSE vectlegende = {minmax:[mini, maxi], unite:varunit}
251
252
253sortie:
254   if keyword_set(key_performance) NE 0 THEN print, 'temps ajoutvect', systime(1)-tempsun
255   return
256end
257
258
Note: See TracBrowser for help on using the repository browser.