; $Id$
;
; Copyright (c) 1983-1998, Research Systems, Inc.  All rights reserved.
;	Unauthorized reproduction prohibited.
; Changes to original routine velovect by Martin Schultz and Franz Rohrer.
;
; Changes to this improved routine by wd (Wolfgang.Dobler@ncl.ac.uk)

PRO WDVELOVECT,U_in,V_in,X_in,Y_in, $
        Missing = Missing, Length = length, Dots = dots,  $
        Color=color, CLIP=clip, NOCLIP=noclip, OVERPLOT=overplot,  $
        NOZERO=nozero, CENTER=center, XRANGE=xrange, YRANGE=yrange,  $
        POSITION=position, POLAR=POLAR, LINEONLY=LINEONLY, $
        ScaLength = scalength, MAXVEC=maxvec, NOPLOT=noplot, QUIET=quiet, $
        _EXTRA=extra
;
;+ 
; NAME:
;	WDVELOVECT
;
; PURPOSE:
;	Produce a two-dimensional velocity field plot.
;
;	A directed arrow is drawn at each point showing the direction and 
;	magnitude of the field.
;
;       This routine was developed from the original VELOVECT routine
;       and now allows two dimensional arrays for X and Y for irregularily
;       spaced data.
;               
; CATEGORY:
;	2D Graphics
;
; CALLING SEQUENCE:
;	WDVELOVECT, U, V [, X, Y]
;
; INPUTS:
;	U:	The X component of the two-dimensional field.  
;		U must be a two-dimensional array.
;
;	V:	The Y component of the two dimensional field.  Y must have
;		the same dimensions as X.  The vector at point [i,j] has a 
;		magnitude of:
;
;			(U[i,j]^2 + V[i,j]^2)^0.5
;
;		and a direction of:
;
;			ATAN2(V[i,j],U[i,j]).
;
;  OR, if the keyword POLAR is set then:
;
;	U:	The magnitude of the two-dimensional field.  
;		U must be a two-dimensional array.
;
;	V:	The angular direction of the two dimensional field (in radians).
;               Y must have same dimensions as X.  
;               Angles are measured as per the computing standard i.e. 
;               anti-clockwise from the +ve x-axis.
;
; OPTIONAL INPUT PARAMETERS:
; 	X:	Optional abcissae values.  X must be a vector with a length 
;		equal to the first dimension of U and V *OR* a 2-dimensional
;               array with the same dimensions as U and V.
;
;	Y:	Optional ordinate values.  Y must be a vector with a length
;		equal to the first dimension of U and V *OR* a 2-dimensional
;               array with the same dimensions as U and V.
;
; KEYWORD INPUT PARAMETERS:
;	COLOR:	The color index used for the plot.
;
;	DOTS:	Set this keyword to 1 to place a dot at each missing point. 
;		Set this keyword to 0 or omit it to draw nothing for missing
;		points.  Has effect only if MISSING is specified.
;
;        POLAR: Use the U and V specified as the magnitude and direction,
;               respectively, of the two-dimensional vector field. 
;               (See U and V above)
;
;
;	LENGTH:	Length factor.  The default of 1.0 makes the longest (U,V)
;		vector the length of a cell.
;
;       SCALENGTH:  Franz Rohrer's modification of the LENGTH keyword:
;              SCALENGTH applies a scale factor relative to the data
;              values, not just a relative scaling.
;              [wd: renamed this from LENGTH. The new behaviour is
;              counterintuitive and certainly should never have
;              changed the behaviour of LENGTH.]
;
;       MISSING: Missing data value.  Vectors with a LENGTH greater
;		than MISSING are ignored.
;
;       OVERPLOT: Set this keyword to make VELOVECT "overplot".  That is, the
;               current graphics screen is not erased, no axes are drawn, and
;               the previously established scaling remains in effect.
;
;       NOZERO: Do not plot zero vectors as dots.
;
;       CENTER: Position of the grid point on the arrow (CENTER=0
;               corresponds to the default behavior of velovect,
;               CENTER=0.5 centers the arrow w.r.t. the grid point and
;               is the default here).
;
;       MAXVEC: Maximum number of arrows; either a scalar (limit total
;               number) or an array [maxvecx,maxvecy]. A value of zero
;               means no limit. The default tries to keep the arrows
;               visible, based on screen resolution.
;
;       NOPLOT: Do not plot any arrows. Useful for just establishing
;               the coordinate system like PLOT_3D_VECT does.
;
;       QUIET:  Don't print informational messages
;
;	Note:   All other keywords are passed directly to the PLOT procedure
;		and may be used to set option such as TITLE, POSITION, 
;		NOERASE, etc.
; OUTPUTS:
;	None.
;
; COMMON BLOCKS:
;	None.
;
; SIDE EFFECTS:
;	Plotting on the selected device is performed.  System
;	variables concerning plotting are changed.
;
; RESTRICTIONS:
;	None.
;
; PROCEDURE:
;	Straightforward.  Unrecognized keywords are passed to the PLOT
;	procedure.  
;
; MODIFICATION HISTORY:
;	DMS, RSI, Oct., 1983.
;	For Sun, DMS, RSI, April, 1989.
;	Added TITLE, Oct, 1990.
;	Added POSITION, NOERASE, COLOR, Feb 91, RES.
;	August, 1993.  Vince Patrick, Adv. Visualization Lab, U. of Maryland,
;		fixed errors in math.
;	August, 1993. DMS, Added _EXTRA keyword inheritance.
;	January, 1994, KDB. Fixed integer math which produced 0 and caused
;		            divide by zero errors.
;	December, 1994, MWR. Added _EXTRA inheritance for PLOTS and OPLOT.
;	June, 1995, MWR. Removed _EXTRA inheritance for PLOTS and changed
;			 OPLOT to PLOTS.
;       September, 1996, GGS. Changed denominator of x_step and y_step vars. 
;       February, 1998, DLD.  Add support for CLIP and NO_CLIP keywords.
;       June, 1998, DLD.  Add support for OVERPLOT keyword.
;       16 Sep 1999: Martin Schultz added support for 2D U and V arrays
;              cleaned up the routine some and added the NOZERO keyword.
;              (renamed as msvelovect.pro on Sep 23)
;              Also included Franz Rohrer's modification of the LENGTH keyword:
;              LENGTH now applies a scale factor relative to the data values,
;              not just a relative scaling. If you don't specify LENGTH it acts
;              as before.
;       29 Feb 2000: mgs. Bug fix: for loop needs long int!
;              Also sx and sy are now always computed even if no x or
;              y are passed into routine.
;       26 Jun 2001: wd. Center arrows and add CENTER keyword
;       27 Jun 2001: wd. Made angles of arrow heads device independent
;        4 Jul 2001: wd. Intercept keywords [XY]RANGE
;       19 Jul 2001: wd. Apply reform() to arguments U and V
;       22 Mar 2002: wd. Reverted F. Rohrer's silly LENGTH behaviour to original
;       22 Mar 2002: wd. Added MAXVEC, QUIET and NOPLOT keywords
;       07 Aug 2004: Antony Mee, U. of Newcastle upon Tyne
;              Added POLAR keyword to allow specification of vectors as
;              magnitude and direction.  Also added LINEONLY to prevent
;              drawing an arrow head.
;-
;
        on_error,2                      ;Return to caller if an error occurs
;
        if (n_elements(quiet) le 0) then quiet = 0
        if (n_elements(noplot) le 0) then noplot = 0
;
        u = reform(u_in)
        v = reform(v_in)
        s = size(u)
        t = size(v)
        if s[0] ne 2 then begin 
baduv:     message, 'U and V parameters must be 2D and same size.'
        endif
        if total(abs(s[0:2]-t[0:2])) ne 0 then goto,baduv
;
        if n_params() lt 3 then begin
           x = findgen(s[1]) 
           sx = size(x)
        endif else begin
           x = x_in
           sx = size(x)
           if (sx[0] eq 2) then begin
              if total(abs(sx[0:2]-s[0:2])) ne 0 then begin
badx:            message, 'X array has incorrect size.'
              endif
           endif else $
              if n_elements(x) ne s[1] then goto,badx
        endelse
;
        if n_params() lt 4 then begin
           y = findgen(s[2]) 
           sy = size(y)
        endif else begin
           y = y_in
           sy = size(y)
           if (sy[0] eq 2) then begin
              if (sx[0] ne 2) then goto,bady
              if total(abs(sy[0:2]-s[0:2])) ne 0 then begin
bady:            message, 'Y array has incorrect size.'
              endif
           endif else $
              if n_elements(y) ne s[2] then goto,bady
        endelse
;
; wd: regrid data if too many points
;
        if (not keyword_set(overplot)) then begin
          ;---------  pretend to plot to get plot window size right  ---------
          if (n_elements(position) lt 4) then begin
            plot, x, y, /NODATA, /NOERASE,XSTYLE=4, YSTYLE=4, TITLE=''
          endif else begin
            plot, x, y, /NODATA, /NOERASE,XSTYLE=4, YSTYLE=4, TITLE='', $
                POSITION=position
          endelse
        endif
        ; only now we can extract position
        position=[!x.window[0],!y.window[0],!x.window[1],!y.window[1]]
        xwidth = !x.window[1] - !x.window[0]
        ywidth = !y.window[1] - !y.window[0]
        case n_elements(maxvec) of
          ;
          0: begin
            ;; default settings: maxvec depends on resolution
            if ((!d.flags and 1) eq 0) then begin ; fixed-size pixels
              npix = 9     ; make longest arrows that many pixels long
              maxvec = [!d.x_vsize*xwidth, !d.y_vsize*ywidth]/npix
            endif else begin    ; scalable pixels
              nmms = 2.0      ;  make longest arrows that many mm long
              maxvec = [!d.x_vsize*xwidth, !d.y_vsize*ywidth]/(nmms*100)
            endelse
          end
          ;
          1: begin
            ;; One arg given : replicate for both directions according
            ;; to number of points in x and y
            maxvec = floor([sqrt(maxvec*sx/sy), sqrt(maxvec*sy/sx)])
          end
          ;
          else:                 ; nothing to do
          ;
        endcase

        regrid = 0
        nx = s[1]
        ny = s[2]
        if ( (maxvec[0] gt 0) and (s[1] gt maxvec[0])) then begin
          nx = maxvec[0]
          regrid = 1
        endif
        if ( (maxvec[1] gt 0) and (s[2] gt maxvec[1])) then begin
          ny = maxvec[1]
          regrid = 1
        endif
        nx = floor(nx)
        ny = floor(ny)
        if (regrid) then begin
          if ((quiet eq 0) and (noplot eq 0)) then begin
            message, /INFO, $
                'Regridding from ' + strtrim(s[1],2) + 'x' + strtrim(s[2],2) $
                + ' to ' + strtrim(nx,2) + 'x' + strtrim(ny,2)
          endif
          x = linspace(minmax(x),nx,GHOST=0.5)
          y = linspace(minmax(y),ny,GHOST=0.5)
          u = congrid(u, nx, ny, /CUBIC)
          v = congrid(v, nx, ny, /CUBIC)
          ; Recalculate these:
          s = size(u)
          t = size(v)
          sx = size(x)
          sy = size(y)
        endif
;
        if (n_elements(center) le 0) then center = 0.5
;
        if n_elements(missing) le 0 then missing = 1.0e30
; ### FR: use LENGTH differently -- allows absolute scaling
;       if n_elements(length) le 0 then length = 1.0

        if keyword_set(POLAR) then begin
          mag = u
        endif else begin
          mag = sqrt(u^2.+v^2.)             ;magnitude.
        endelse
                ;Subscripts of good elements
        nbad = 0                        ;# of missing points
; ## mgs: because of defaulting 5 lines above, missing always has a value!!
;       if n_elements(missing) gt 0 then begin
           good = where(mag lt missing) 
           if keyword_set(dots) then bad = where(mag ge missing, nbad)
;       endif else begin
;               good = lindgen(n_elements(mag))
;       endelse
        ugood = u[good]
        vgood = v[good]
        x0 = min(x)                     ;get scaling
        x1 = max(x)
        y0 = min(y)
        y1 = max(y)
	x_step=(x1-x0)/(s[1]-1.0)   ; Convert to float. Integer math
	y_step=(y1-y0)/(s[2]-1.0)   ; could result in divide by 0

        if keyword_set(POLAR) then begin
	  maxmag=max(abs(u))/min([x_step,y_step])
;AJWM:   Strange but true!
;        sina and cosa are switched to use angles measured CCW from the +ve 
;        +ve x-axis switching them back gives angles measured CW from 
;        the +ve y-axis; the (x,y)-component based code considers the
;        angles of vectors such.
;
;        Could consider some keywords to allow switching between the two
;        methods.  But will stick with the IDL/computing standard for now.
;
          default,length,1.0
          cosa = length*ugood/maxmag*sin(vgood)
          sina = length*ugood/maxmag*cos(vgood)
        endif else begin
	  maxmag=max([max(abs(ugood/x_step)),max(abs(vgood/y_step))])
; ### FR:
          if n_elements(scalength) gt 0 then maxmag=scalength/x_step
          sina = (ugood/maxmag)
          cosa = (vgood/maxmag)
; ### wd: LENGTH overrides SCALENGTH
          if (n_elements(length) gt 0) then begin
            sina = length * (ugood/maxmag)
            cosa = length * (vgood/maxmag)
          endif
        endelse
;
        if n_elements(title) le 0 then title = ''
        ;--------------  plot to get axes  ---------------
        if n_elements(color) eq 0 then color = !p.color
        if n_elements(noclip) eq 0 then noclip = 0
        x_b0=x0-x_step
	x_b1=x1+x_step
	y_b0=y0-y_step
	y_b1=y1+y_step
        if (not keyword_set(overplot)) then begin
;         if n_elements(position) eq 0 then begin
          ;; The following would be clearer in the form
          ;;   if ((n_elts(xr) le 1) or (xr[0] eq xr[1])) then xr = [..] ,
          ;; but IDL doesn't support logical shortcircuiting
          if (n_elements(xrange) le 1) then xrange=[x_b0,x_b1]
          if (xrange[0] eq xrange[1])  then xrange=[x_b0,x_b1]
          if (n_elements(yrange) le 1) then yrange=[y_b0,y_b1]
          if (yrange[0] eq yrange[1])  then yrange=[y_b0,y_b1]
          plot,[x_b0,x_b1],[y_b1,y_b0],/nodata,/xst,/yst, $
              color=color, xrange=xrange, yrange=yrange, POSITION=position, $
              _EXTRA = extra
;         endif else begin
;           plot,[x_b0,x_b1],[y_b1,y_b0],/nodata,/xst,/yst, $
;             color=color, _EXTRA = extra
;         endelse
        endif
; Only plot if NOPLOT is not set
        if (noplot eq 0) then begin
            if n_elements(clip) eq 0 then $
                clip = [!x.crange[0],!y.crange[0],!x.crange[1],!y.crange[1]]
;
            r = .3                ;len of arrow head
            angle = 22.5 * !dtor  ;Angle of arrowhead
            st = r * sin(angle)   ;sin 22.5 degs * length of head
            ct = r * cos(angle)
;
; Make angles of arrow head device independent:
            xd=(!x.crange[1]-!x.crange[0])/(!x.window[1]-!x.window[0])/!d.x_size
            yd=(!y.crange[1]-!y.crange[0])/(!y.window[1]-!y.window[0])/!d.y_size
;
            for i=0L,n_elements(good)-1 do begin ;Each point
                if (sx[0] eq 2) then begin
                    x0 = x[good[i]] ;get coords of start & end
                    y0 = y[good[i]]
                endif else begin
                    x0 = x[good[i] mod s[1]] ;get coords of start & end
                    y0 = y[good[i] / s[1]]
                endelse
                dx = sina[i]
                x0 = x0 - center*dx
                x1 = x0 + dx
                dy = cosa[i]
                y0 = y0 - center*dy
                y1 = y0 + dy
;	         xd=x_step
;	         yd=y_step
                ; plot zero vectors as dots
                if (mag[good[i]] eq 0.) then begin
                    if  (not keyword_set(NOZERO)) then $
                        plots,x0,y0,psym=3,color=color,clip=clip,  $
                        noclip=noclip  
                endif else $
                  if keyword_set(LINEONLY) then begin
                    plots,[x0,x1],[y0,y1], $
                           color=color,clip=clip,noclip=noclip
                  endif else begin
                    plots,[x0,x1,x1-(ct*dx/xd-st*dy/yd)*xd, $
                           x1,x1-(ct*dx/xd+st*dy/yd)*xd], $
                          [y0,y1,y1-(ct*dy/yd+st*dx/xd)*yd, $
                           y1,y1-(ct*dy/yd-st*dx/xd)*yd], $
                           color=color,clip=clip,noclip=noclip
                  endelse
            endfor
            if nbad gt 0 then begin
                if (sx[0] eq 2) then begin ;Dots for missing?
                    PLOTS, x[bad], y[bad], psym=3, $
                        color=color, clip=clip,noclip=noclip   
                endif else begin
                    PLOTS, x[bad mod s[1]], y[bad / s[1]], psym=3, $
                        color=color, clip=clip,noclip=noclip
                endelse
            endif
        endif
;
end