PRO plot_vel_image, array, velmax=velmax,_extra=_extra, white_val=white_val, $
                    black_val=black_val, reverse=reverse

;+
; NAME
;
;     PLOT_VEL_IMAGE
;
; EXPLANATION
;
;     A version of PLOT_IMAGE for velocity maps that implicitly sets
;     the color table to display reds and blues. The velocity image
;     has white as the zero velocity point, and also the axes are
;     colored white or black as keeping with normal plot_image
;     behavior. This is achieved by setting a color table such that:
;     bit 0 is reserved for black; bits 254 and 255 are white; bits
;     1-126 are for blueshifts; bit 127 is zero shift (white) and bits
;     128-253 are for redshifts.
;
;     This routine will accept any of the keywords accepted by
;     PLOT_IMAGE. 
;
; INPUTS
;
;     ARRAY   A 2-D array to be plotted.
;
; OPTIONAL INPUTS
;
;     VELMAX  The image will be scaled so that values will lie between
;             +/- VELMAX. Default value is 30 km/s.
;
;     WHITE_VAL Any pixels that have a value of WHITE_VAL will be
;               colored white in the plotted image. This is intended
;               for distinguishing missing data points in the velocity
;               map. 
;
;     BLACK_VAL Any pixels that have a value of BLACK_VAL will be
;               colored black in the plotted image. This is intended
;               for distinguishing missing data points in the velocity
;               map. 
;
; HISTORY
;
;     Ver.1, 6-May-2011, Peter Young
;-

IF n_params() LT 1 THEN BEGIN
  print,'Use:  IDL> plot_vel_image, array, velmax=velmax'
  return
ENDIF 

IF n_elements(velmax) EQ 0 THEN velmax=30

n_w=0
IF n_elements(white_val) NE 0 THEN BEGIN
  i_w=where(array EQ white_val,n_w)
ENDIF

n_b=0
IF n_elements(black_val) NE 0 THEN BEGIN
  i_b=where(array EQ black_val,n_b)
ENDIF 

arr=byte(findgen(126)*253/125)

r=bytarr(256)
r[0]=0
r[127]=255
r[128:253]=255
r[1:126]=arr
r[254:255]=255
;
g=bytarr(256)
g[0]=0
g[127]=255
g[1:126]=arr
g[128:253]=reverse(arr)
g[254:255]=255
;
b=bytarr(256)
b[0]=0
b[127]=255
b[128:253]=reverse(arr)
b[1:126]=255
b[254:255]=255

IF keyword_set(reverse) THEN BEGIN
  rx=r
  r=b
  b=rx
ENDIF 

tvlct,r,g,b

array_byte=bytscl(array,min=-velmax,max=velmax,top=252)
array_byte=array_byte+1b

IF n_b NE 0 THEN array_byte[i_b]=0
IF n_w NE 0 THEN array_byte[i_w]=255

plot_image,array_byte,/noscale,_extra=_extra

END