[2] | 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 | ;------------------------------------------------------------ |
---|
| 55 | pro 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 | ;----------------------------------------------------------- |
---|
| 65 | ;------------------------------------------------------------ |
---|
| 66 | ; on recupere les eventuelles info sur les champs |
---|
| 67 | ;------------------------------------------------------------ |
---|
| 68 | grilleu = litchamp(vecteur.(0), /grid) |
---|
| 69 | if grilleu EQ '' then grilleu = 'U' |
---|
| 70 | grillev = litchamp(vecteur.(1), /grid) |
---|
| 71 | if grillev EQ '' then grillev = 'V' |
---|
[41] | 72 | |
---|
[2] | 73 | IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1 |
---|
| 74 | IF grilleu EQ grillev THEN interpolle = 0 ELSE interpolle = 1 |
---|
| 75 | if keyword_set(inverse) then begin |
---|
| 76 | rien = u |
---|
| 77 | u = v |
---|
| 78 | v = rien |
---|
| 79 | endif |
---|
| 80 | ;------------------------------------------------------------ |
---|
| 81 | ; on trouve les points que u et v ont en communs |
---|
| 82 | ;------------------------------------------------------------ |
---|
| 83 | if interpolle then begin |
---|
[41] | 84 | indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] |
---|
| 85 | indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] |
---|
[2] | 86 | indicex = inter(indicexu, indicexv) |
---|
[41] | 87 | indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] |
---|
| 88 | indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] |
---|
[2] | 89 | indicey = inter(indiceyu, indiceyv) |
---|
| 90 | nx = n_elements(indicex) |
---|
| 91 | ny = n_elements(indicey) |
---|
| 92 | indice2d = lindgen(jpi, jpj) |
---|
| 93 | indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] |
---|
| 94 | ;------------------------------------------------------------ |
---|
| 95 | ; extraction de u et v sur le domaine qui convient |
---|
| 96 | ;------------------------------------------------------------ |
---|
| 97 | case 1 of |
---|
| 98 | (size(u))[0] NE 2 OR (size(v))[0] NE 2: return |
---|
| 99 | (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ |
---|
| 100 | (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN |
---|
| 101 | if nxu NE nx then $ |
---|
[41] | 102 | if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] |
---|
[2] | 103 | IF nxv NE nx THEN $ |
---|
[41] | 104 | if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] |
---|
[2] | 105 | IF nyu NE ny THEN $ |
---|
[41] | 106 | if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] |
---|
[2] | 107 | IF nyv NE ny THEN $ |
---|
[41] | 108 | if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] |
---|
[2] | 109 | END |
---|
| 110 | (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ |
---|
| 111 | (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN |
---|
| 112 | u = u[indice2d] |
---|
| 113 | v = v[indice2d] |
---|
| 114 | END |
---|
| 115 | ELSE:BEGIN |
---|
| 116 | ras = report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') |
---|
| 117 | return |
---|
| 118 | end |
---|
| 119 | endcase |
---|
| 120 | ;------------------------------------------------------------------ |
---|
| 121 | ; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee |
---|
| 122 | ;------------------------------------------------------------------ |
---|
| 123 | if ny EQ 1 then begin |
---|
| 124 | u = reform(u, nx, ny) |
---|
| 125 | v = reform(v, nx, ny) |
---|
| 126 | endif |
---|
| 127 | ;------------------------------------------------------------------ |
---|
| 128 | ; construction de u et v aux pts T |
---|
| 129 | ;----------------------------------------------------------- |
---|
| 130 | a=u(0,*) |
---|
| 131 | u=(u+shift(u,1,0))/2. |
---|
[41] | 132 | if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*)=a |
---|
[2] | 133 | a=v(*,0) |
---|
| 134 | v=(v+shift(v,0,1))/2. |
---|
[41] | 135 | if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0)=a |
---|
[2] | 136 | ;---------------------------------------------------------------------------- |
---|
| 137 | ; attribution du mask et des tableau de longitude et latitude |
---|
| 138 | ; on recupere la grille complette pour etablir un grand mask etendu ds les 4 |
---|
| 139 | ; directions pour couvrir les points pour lesquels un pt terre a ete pris en |
---|
| 140 | ; compte (faire un petit dessin...) |
---|
| 141 | ;---------------------------------------------------------------------------- |
---|
| 142 | vargrid='T' |
---|
[41] | 143 | msku = (umask())[indice2d+jpi*jpj*firstzt] |
---|
| 144 | mskv = (vmask())[indice2d+jpi*jpj*firstzt] |
---|
[2] | 145 | glam = glamt[indice2d] |
---|
| 146 | gphi = gphit[indice2d] |
---|
| 147 | if ny EQ 1 then begin |
---|
| 148 | msku = reform(msku, nx, ny) |
---|
| 149 | mskv = reform(mskv, nx, ny) |
---|
| 150 | ; glam = reform(glam, nx, ny) |
---|
| 151 | ; gphi = reform(gphi, nx, ny) |
---|
| 152 | endif |
---|
| 153 | ;----------------------------------------------------------- |
---|
| 154 | ; on masque u et v le long des cotes (la on l''on ne peut pas calculer |
---|
| 155 | ; la moyenne) |
---|
| 156 | ;----------------------------------------------------------- |
---|
| 157 | ; extention du mask |
---|
| 158 | u = u*msku*shift(msku,1,0) |
---|
| 159 | v = v*mskv*shift(mskv,0,1) |
---|
| 160 | ENDIF ELSE BEGIN |
---|
[41] | 161 | u = u*tmask[firstxt:lastxt,firstyt:lastyt,firstzt] |
---|
| 162 | v = v*tmask[firstxt:lastxt,firstyt:lastyt,firstzt] |
---|
[2] | 163 | indice2d = lindgen(jpi, jpj) |
---|
[41] | 164 | indice2d = indice2d[firstxt:lastxt, firstyt:lastyt] |
---|
[2] | 165 | nx = nxt |
---|
| 166 | ny = nyt |
---|
| 167 | endelse |
---|
| 168 | tabnorme=sqrt(u^2+v^2) |
---|
| 169 | nan = where(finite(u, /nan) EQ 1) |
---|
| 170 | if nan[0] NE -1 then u[nan] = 1e5 |
---|
| 171 | nan = where(finite(v, /nan) EQ 1) |
---|
| 172 | if nan[0] NE -1 then v[nan] = 1e5 |
---|
| 173 | if keyword_set(vectmin) then BEGIN |
---|
[41] | 174 | toosmall=where(tabnorme lt vectmin) |
---|
| 175 | if toosmall[0] NE -1 then begin |
---|
| 176 | u[toosmall] = 1e5 |
---|
| 177 | v[toosmall] = 1e5 |
---|
[2] | 178 | ENDIF |
---|
| 179 | endif |
---|
| 180 | if keyword_set(vectmax) then BEGIN |
---|
[41] | 181 | toobig=where(tabnorme gt vectmax) |
---|
| 182 | if toobig[0] NE -1 then begin |
---|
| 183 | u[toobig] = 1e5 |
---|
| 184 | v[toobig] = 1e5 |
---|
[2] | 185 | ENDIF |
---|
| 186 | ENDIF |
---|
| 187 | ;----------------------------------------------------------- |
---|
| 188 | ; remise d''une grande valeur sur tous les points pour lesquelles on ne |
---|
| 189 | ; put faire le calcul |
---|
| 190 | ;----------------------------------------------------------- |
---|
| 191 | if interpolle then t2 = msku*shift(msku,1,0)*mskv*shift(mskv,0,1) $ |
---|
[41] | 192 | ELSE t2 = tmask[firstxt:lastxt,firstyt:lastyt,firstzt] |
---|
| 193 | if NOT keyword_set(key_periodic) OR nx NE jpi then t2(0, *)=0. |
---|
[2] | 194 | t2(*,0)=0. |
---|
| 195 | terre=where(t2 eq 0) |
---|
| 196 | if terre[0] ne -1 then begin |
---|
| 197 | u(terre)=1e5 |
---|
| 198 | v(terre)=1e5 |
---|
| 199 | ENDIF |
---|
| 200 | ;----------------------------------------------------------- |
---|
| 201 | ; tracer qu'un vecteur sur |
---|
| 202 | ;----------------------------------------------------------- |
---|
| 203 | if keyword_set(unvectsur) then BEGIN ; |
---|
| 204 | ; indx est un vecteur contenant les numero des colonnes a selectionner |
---|
| 205 | ; indy est un vecteur contenant les numero des lignes a selectionner |
---|
| 206 | if n_elements(unvectsur) EQ 1 then begin |
---|
| 207 | indx = where((lindgen(nx) MOD unvectsur[0]) eq 0) |
---|
| 208 | indy = where((lindgen(ny) MOD unvectsur[0]) eq 0) |
---|
| 209 | ENDIF ELSE BEGIN |
---|
| 210 | indx = where((lindgen(nx) MOD unvectsur[0]) eq 0) |
---|
| 211 | indy = where((lindgen(ny) MOD unvectsur[1]) eq 0) |
---|
| 212 | ENDELSE |
---|
| 213 | ; a partir de indx et indy on va construire un tableau d''indices 2d |
---|
| 214 | ; qui donnera les indices des points intersections des colonnes |
---|
| 215 | ; specifiee par indx |
---|
| 216 | indicereduit = indx#replicate(1,n_elements(indy))+nx*replicate(1,n_elements(indx))#indy |
---|
| 217 | ; on reduit les tableaux qui vont etre passes a vecteur. |
---|
| 218 | u = u[indicereduit] |
---|
| 219 | v = v[indicereduit] |
---|
| 220 | tabnorme = tabnorme[indicereduit] |
---|
| 221 | ; |
---|
| 222 | endif |
---|
| 223 | ;----------------------------------------------------------- |
---|
| 224 | ; |
---|
| 225 | ;----------------------------------------------------------- |
---|
| 226 | if keyword_set(inverse) then begin |
---|
| 227 | rien = u |
---|
| 228 | u = v |
---|
| 229 | v = rien |
---|
| 230 | endif |
---|
| 231 | ;----------------------------------------------------------- |
---|
| 232 | ; trace des vecteurs |
---|
| 233 | ;---------------------------------------------------------- |
---|
| 234 | vecteur, u, v, tabnorme, indice2d, indicereduit, missing=1e5, _extra = ex |
---|
| 235 | ;----------------------------------------------------------- |
---|
| 236 | ; on complete la legende |
---|
| 237 | ;----------------------------------------------------------- |
---|
| 238 | if terre[0] ne -1 then mini = min(tabnorme[where(t2 eq 1)], max = maxi, /nan) $ |
---|
| 239 | ELSE mini = min(tabnorme, max = maxi, /nan) |
---|
| 240 | |
---|
| 241 | if litchamp(vecteur.(0), /u) NE '' then $ |
---|
| 242 | vectlegende = {minmax:[mini, maxi], unite:litchamp(vecteur.(0), /u)} $ |
---|
| 243 | ELSE vectlegende = {minmax:[mini, maxi], unite:varunit} |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | sortie: |
---|
| 247 | if keyword_set(key_performance) NE 0 THEN print, 'temps ajoutvect', systime(1)-tempsun |
---|
| 248 | return |
---|
| 249 | end |
---|
| 250 | |
---|
| 251 | |
---|