;+------------------------------------------------------------
; NAME:
;       TH_IMAGE_CONT
; PURPOSE:
;       Plot contours of Images with both contour lines and colors
; CATEGORY:
;	2-D Plotting
; USAGE:
;******************** HIGH FREQUENCY RADAR DIVISION, SRL **********************
;*************************** Ionospheric Effects ******************************
;HELP
;1 TH_IMAGE_CONT
;	Overlays an image and a contour plot and optionally adds a scale bar.
;	Based on the IDL USERLIB routine IMAGE_CONT. This routine supersedes 
;	the USERLIB one, having far more functionality, yet capable of EXACTLY
;	reproducing the effect of IMAGE_CONT. Scale bar appears on the 
;	right-hand-side unless /NOBAR is set. NB: the scale bar is
;	automatically created by a recursive call to this routine using the
;	same colour and image parameters 
;
; Format:
;	In its simplest form (allowing all parameters to default)
;
;    IDL> TH_IMAGE_CONT, IMAGE
;
;	or,
;
;    IDL> TH_IMAGE_CONT, IMAGE, X, Y
;
;	And in its most complex form, specifying ALL parameter
;
;    IDL> TH_IMAGE_CONT, IMAGE, $
;		ASPECT=aspect, $
;		BADPTS=badpts,$
;     		BAR_LABEL_LENGTH=bar_label_length, $
;		BAR_RANGE=bar_range, $
;		BAR_SEPARATION=bar_separation, $
;		BAR_TICKLEN=bar_ticklen, $
;		BAR_TICKNAME=bar_tickname, $
;		BAR_TICKS=bar_ticks, $
;		BAR_TICKV=bar_tickv,$
;		BARSZ_CHARS=barsz_chars, $
;		BOTTOMCOLOUR=bottomcolour,$
;		C_COLORS=c_colors,$
;		C_LINESTYLE=c_linestyle, $
;		C_THICK=c_thick, $
;		CONGRID=congrid, $
;		CONT=cont, $
;		CRANGE=crange, $
;		CT=ct, $
;		CUBIC=cubic, $
;		DEBUG=debug, $
;		EXACT=exact, $
;		RAISE_PTITLE=raise_ptitle
;		IMAGE_WINDOW=image_window, $
;		INTERP=interp,$
;		LEVELS=levels,$
;		MAX_VALUE=max_value, $
;		NLEVELS=nlevels, $
;		NOBAR=nobar, $
;		NOCONT=nocont, $
;		NOERASE=noerase,$
;		TOPCOLOUR=topcolour,$
;		TSIZE=tsize, $
;		WINDOW_SCALE=window_scale, $
;		XRANGE=xrange, $
;		YRANGE=yrange
; 
;2	IMAGE
;		2-dimensional array to display as an image.
;
;2	/ASPECT
;		set to retain image's aspect ratio.  Assumes square
;		pixels.  If /WINDOW_SCALE is set, the aspect ratio is retained.
;
;2	BADPTS
;		indices into IMAGE data which define the bad points.
;		These will not be contoured
;
;2 BAR_LABEL_LENGTH
;     the space, in characters (default = 2), for the bar label
;2	BAR_RANGE	
;	 	set the range limits for the colour scale bar
;		(same as CRANGE, defaults to autoscaling if
;		 BAR_TICKV not set)
;
;2	BAR_SEPARATION	
;		the separation, in characters (default = 2), between the scale
;		bar and the image. Note that the |y-ticklength| will be added
;		to this value if y-ticklength < 0.
;		Both the image and the colour bar need to fit 
;		into the space allowed for the plot window, otherwise an
;		informative message will be printed and unpredictable results
;		may occur for the displayed image 
;
;2	BAR_TICKS
;		set the number of tick intervals for the labelling
;		of the scale bar 
;		(defaults to !z.ticks)
;
;2	BAR_TICKV
;		the values to label on the scale bar. If this is set
;		then the scale bar will have AT LEAST this range
;		(defaults to !z.tickv)
;
;2	BAR_TICKNAME
;		the labels to use on the scale bar 
;		(defaults to !z.range)
;
;2	BAR_TICKLEN
;		the length of the ticks on the scale bar  in fractions of tick
;		bar window, (defaults to !z.ticklen)
;
;2	BARSZ_CHARS
;		the size of the scale bar in characters (default = 2). 
;		Both the image and the colour bar need to fit 
;		into the space allowed for the plot window, otherwise an
;		informative message will be printed and unpredictable results
;		may occur for the displayed image. If the value of this
;		keyword is <=0 then no bar will be displayed BUT the scaling
;		and window size will be calculated as those a colour bar is to
;		be used. This is useful when doing multiple plots per page
;		where only some scale bars are not required but you want the
;		plots to all be the same size. Set BARSZ_CHARS = -#chars to
;		allow room for a bar of #chars size but not to put a scale bar
;		on the plot. Then set BARSZ_CHARS = +#chars to plot the bar on
;		alongside another plot, to end up with images of the same
;		size. Useful in collaboration with the SIDES procedure (which
;		will set flags for when the plot is on the Left,Right,Top and
;		Bottom of the plot window).
;
;
;2	BOTTOMCOLOUR
;		Set this keyword to the colour index of the desired bottom
;		colour (range from 0 to TOPCOLOUR-1).
;		Note that the default value for this keyword is 1,
;		which allows the colour of the image to
;		be independent of the background and axes colours
;		(!P.background and !P.color). If the user sets this keyword
;		then allowance should be made for these colours as they are
;		generally swapped for POSTSCRIPT and X-Window devices
;
;2	/CONGRID 	
;		if the image has to be resampled then use the USERLIB CONGRID
;		routine 
;
;2	/CONT 	
;		only do the contouring (no image)
;
;2	CRANGE	
;		set the range limits for the colour scale bar
;		(same as BAR_RANGE, defaults to autoscaling if
;		 BAR_TICKV not set)
;
;2	CT	
;		load a color table (uses LOADCT via MK_COLOR)
;
;2	/CUBIC	
;		if the image has to be resampled AND interpolated then use the
;		CUBIC interpolation rather than the bi-linear (see INTERP
;		Keyword) 
;
;2	/DEBUG	
;		write out some inforamtion as it goes
;
;
;2	/EXACT
; 
;		When set this will force the contour routine to fit to the
;		exact positions relative to the image.
;		When data is displayed as an image each datum is expanded out
;		to fill a pixel of finite dimensions. The assignment of where
;		the "data value" resides within this space is open to debate,
;		but is most appropriately (to this author) assigned to the
;		geometric centre of the pixel. Most defaults assign this
;		position to be at the bottom left-hand corner of the
;		pixel. Contour will fit to the 2-d plane assuming that the
;		data value is associated with the bottom left-hand corner. To
;		reconcile this with the notion of the value being in the
;		middle of the pixel the contour call is made with the x/y
;		values and ranges for the image adjusted (effectively by half
;		a pixel in the x/y directions). This is the EXACT mapping. By
;		default, the mapping will be the default contour one.
;
;2	RAISE_PTITLE
;		Raise the plot title by this many character units above the
;		plot to allow room to put a label on the top x-axis. Default
;		is raise by 1 char. If not called then the default y-position
;		is 1 char unit above plot (allowing room for xticks, and
;		scaled by !P.charsize)
;
;2	IMAGE_WINDOW
;		the position of the image window. Can be used to set
;		!p.position so you can over-plot the image.
;		(Only useful when the scale-bar has been used)
;
;2	/INTERP 
;		set to bi-linear interpolate if image is resampled.
;		(see also the CUBIC keyword)
;
;2	/NOBAR 
;		dont put a scale bar on the right-hand-side
;
;2	/NOCONT
;		only do the imaging (no contours)
;
;2	/NOERASE
;		dont erase the previous plot
;       /NOLOAD
;		if set, do not load color table
;
;2	TOPCOLOUR
;		Set this keyword to the colour index of the desired bottom
;		colour (range from BOTTOMCOLOUR+1 to !D.n_colors-1). Note that
;		the default value for this keyword is !D.n_colors-2, which
;		allows the colour of the image to be independent of the
;		background and axes colours (!P.background and !P.color). If
;		the user sets this keyword then allowance should be made for
;		these colours as they are generally swapped for POSTSCRIPT and
;		X-Window devices 
;
;2	TSIZE	
;		size of the plot title (default = 1)
; 
;2	/WINDOW_SCALE
;		set to scale the window size to the image size,
;		otherwise the image size is scaled to the window size.
;		Ignored when outputting to devices with scalable pixels.
;
;2	XRANGE
;		will set the ranges for the x-axes labelling
;
;2	YRANGE
;		will set the ranges for the y-axes labelling
;
;2	Contour
;		most of the CONTOUR parameters are passed directly
;
;2 Examples
; IDL>   th_image_cont, image
;
; IDL>	th_image_cont, image, /nocont, /nobar
;
; IDL>	!p.title = "!17 This is an example of what can be done"
; IDL>	!x.title = "X Title"	& !y.title = "Y Title"	& !z.title = "Z Title"
; IDL>	!x.ticklen = -0.02	& !y.ticklen = -0.02	& !z.ticklen = -0.02
; IDL>	!p.charsize = 1.5
; IDL>	levels = findgen(5)*2
; IDL>	image = findgen(20,20)/40.
; IDL>	th_image_cont, image, crange=[0,10], /follow, level=levels, $
;			tsize=1.5*!p.charsize, bar_tickv=levels, c_char=1.5
;
;2 Error_responses
; Returns to the calling procedure on an error
;
;2 Limitations/Assumptions 
;	The currently selected display is affected.
;	If the device has scalable pixels then the image is written over
;	the plot window.
;	As with all TV style image displays, the axes range is independent of 
;	the image, so it is up to the user the ensure correct labelling of the
;	axes.
;	NOTE: if the user aborts while this routine is processing then the 
;	system variables (in particular !p.position) will have 
;	been changed, causing subsequent plots to appear different. Issue the 
;	command "resetplt,/all" to reset all the system variables back to the 
;	startup state.
;
;	NOTE: if x or y does not have a constant interval, the image will be
;		wrong. Instead, use contour, /fill...
;2 References 
; 	See USERLIB IMAGE_CONT
;
;2 Keywords
;	Graphics images contours
;
;2 MODIFICATIONS:
;  -------------
;  06-Jun-2006 several changes including adding x & ytickformat keywords
;  06-Apr-2006 limit max dimension of interpolated postscript plots to 2000 pixels
;  14-Mar-2006 added zrange keyword [BD]
;  25-Oct-01 Made Contours work when /ASPECT set. Made Contours work when
;	     XRANGE or YRANGE set.
;  18-Apr-00 Made interpolated image align with data. /EXACT does not give
;		  an exact representation. Warn if style=2
;  06-Dec-99 Added another pixel to TV image size, so it fills the grid square
;  23-Aug-99 Allow Contour, Z, X, Y format. Use MK_COLOR
;  07-Jan-98 BD allow more space at the bottom and right-hand side for text
;               (keyword BAR_LABEL_LENGTH added)
;  07-Aug-97 Bill Davis- v. 2.24
;                [1] added _extra call to main CONTOUR and removed
;                    it from those making the bar.
; (more modifications in file)
;-
;2 Amendments;	
; 	May, 1988, DMS		1.00
;...
;-----------------------------------------------------------------
;  get a subimage from specified x & y or !x.range & !y.range
function subimg, im, x_opt, y_opt

   ; make a subimage, if !x.range & !y.range less than full

IF (!X.range[0] EQ !X.range[1] ) AND $
   (!Y.range[0] EQ !Y.range[1] ) THEN RETURN, im

array_info = SIZE( im )
xsize = array_info( 1 )
ysize = array_info( 2 )
if xsize le 2 or ysize le 2 then return, im

if n_elements( x_opt ) GT 0 then x = x_opt $
                            else x = findgen( xsize )
if n_elements( y_opt ) GT 0 then y = y_opt $
                            else y = findgen( ysize )
ix1 = NEARESTI( x, !X.range[0])
ix2 = NEARESTI( x, !X.range[1])
iy1 = NEARESTI( y, !y.range[0])
iy2 = NEARESTI( y, !y.range[1])

return, im[ ix1>0 : ix2>(ix1+1), iy1>0 : iy2>(iy1+1) ]
end
;**************************** TH_IMAGE_CONT Code *****************************
;-----------------------------------------------------------------
PRO Th_image_cont, image_in, xvalues_opt_in, yvalues_opt_in, $
                   WINDOW_SCALE=window_scale, $
                   ASPECT=aspect, $
                   DEBUG=debug,  $
                   INTERP=interp, $
                   CONT=cont, $
                   NOCONT=nocont,  $
                   IMAGE_WINDOW=image_window, $
                   NOERASE=noerase, $
                   NOBAR=nobar,  $
                   CT=ct, $
                   MAX_VALUE=max_value, $
                   LEVELS=levels,  $
                   C_COLORS=c_colors, $
                   C_LINESTYLE=c_linestyle,  $
                   C_THICK=c_thick, $
                   NLEVELS=nlevels, $
                   XRANGE=xrange,  $
                   YRANGE=yrange, $
                   ZRANGE=zrange, $
                   BADPTS=badpts, $
                   BAR_LABEL_LENGTH=bar_label_length, $
                   BAR_TICKS=bar_ticks, $
                   BAR_TICKNAME=bar_tickname, $
                   BAR_TICKV=bar_tickv, $
                   BAR_TICKLEN=bar_ticklen, $
                   BARSZ_CHARS=barsz_chars, $
                   TSIZE=tsize, $
                   BAR_RANGE=bar_range, $
                   CRANGE=crange,  $
                   CUBIC=cubic, $
                   CONGRID=congrid_keyword, $
                   TOPCOLOUR=topc, $
                   BOTTOMCOLOUR=bottomc, $
                   BAR_SEPARATION=bar_separation, $
                   EXACT=exact, $
                   RAISE_PTITLE=raise_ptitle, $
		   MAX_COLORS_TO_USE=MAX_COLORS_TO_USE, $
		   noLoad=noLoad, $
	  btitle=btitle, bcharsize=bcharsize, $
	  xtickformat=xtickformat, ytickformat=ytickformat,   $
	  currentTime=currentTime, $
	  outPos=outPos, $
                   _EXTRA=_extra
  COMMON Pset, n, lfile, encaps, colour, bppix

  IF (KEYWORD_SET(debug) ) THEN  $
    message, /INFO, ' Entered Procedure......... please wait '

  IF (KEYWORD_SET(debug) ) THEN debug = 1 else debug = 0
  IF debug then BEGIN
     print,'!x.range=', !x.range
;;;     print,'xrange=', xrange
     print,'!y.range=', !y.range
  ENDIF
  
  if not debug then  on_error, 2                   ;Return to caller if an error occurs

  IF (!D.name EQ 'NULL') THEN RETURN ;no more to do
  
  if n_elements( noLoad ) eq 0 then noLoad = 0
  
  if NOT  noLoad then begin
     if n_elements( ct ) eq 0 then ct=3	; hot metal
     if n_elements( MAX_COLORS_TO_USE ) eq 0 then MAX_COLORS_TO_USE = 128

     dum=MK_COLOR(table=ct, N_NONFIXED = n_nonfixed, MAX_COLORS_TO_USE=MAX_COLORS_TO_USE)
  endif else begin
;;;     dum=MK_COLOR( N_NONFIXED = n_nonfixed )
        ; these probably don't matter when /noLoad set:
     if n_elements( MAX_COLORS_TO_USE ) eq 0 then MAX_COLORS_TO_USE = 256
     if n_elements( N_NONFIXED ) eq 0 then N_NONFIXED=256 
  endelse

  orig_p = !P
  orig_x = !X
  orig_y = !Y	
  orig_z = !Z
  IF debug then  print,'!y.range=',!y.range

  image=reform(image_in)

  if n_elements( zrange ) eq 2 then image = zrange[0] > image < zrange[1]
  if n_elements( btitle ) gt 0 then !z.title = btitle
  
  ainfo = size( image )
  
  IF ainfo[0] LT 2 THEN BEGIN
    ;;; message, '1st parameter must be 2D array'
    print, '1st parameter must be 2D array'
     return
  endif
  if n_elements( xvalues_opt_in ) eq 0 then xvalues_opt_in=findgen( ainfo[1] )
  if n_elements( yvalues_opt_in ) eq 0 then yvalues_opt_in=findgen( ainfo[2] )
  xvalues_opt= xvalues_opt_in
  yvalues_opt= yvalues_opt_in

  dx = xvalues_opt[1:*] - xvalues_opt[0:*] 
  negs = where( dx lt 0, nnegs)
  if nnegs gt 0 and nnegs lt n_elements( dx ) then begin
     print, '*** X must be monotonic'
     return
  endif

  dy = yvalues_opt[1:*] - yvalues_opt[0:*] 
  negs = where( dy lt 0, nnegs)
  if nnegs gt 0 and nnegs lt n_elements( dy ) then begin
     print, '*** Y must be monotonic'
     return
  endif

  IF N_Params(0) ge 3 THEN BEGIN
     Use_Args = 1
     if (!X.range[0] eq !X.range[1] ) then  $
         !X.range = [min(xvalues_opt),max(xvalues_opt)]
     if (!Y.range[0] eq !Y.range[1] ) then  $
         !Y.range = [min(yvalues_opt),max(yvalues_opt)]
  ENDIF ELSE Use_Args = 0
  
  IF (KEYWORD_SET(xrange) ) THEN BEGIN
     if n_elements(xrange) eq 1 then xrange=[0,0]
     !X.range = float(xrange)
     congrid_keyword = 1	; can't get poly_2d in this case
  ENDIF
  IF (KEYWORD_SET(yrange) ) THEN  BEGIN
     !Y.range = float(yrange) 
     congrid_keyword = 1	; can't get poly_2d in this case
  ENDIF
  IF (KEYWORD_SET(crange) ) THEN bar_range = crange
  if n_elements( bar_range ) lt 2 then begin
     if n_elements( zrange ) eq 2 then bar_range = zrange
  endif
  
  IF debug then BEGIN
     print,'POS 1.1 !x.range=', !x.range
;;;     print,'POS 1.1 xrange=', xrange
  ENDIF
  
  if Use_args then begin
     x4lim = xvalues_opt
     y4lim = yvalues_opt
  endif else begin
     x4lim = indgen( aInfo[1] )
     y4lim = indgen( aInfo[2] )
  endelse
  
  xinds = where( x4lim ge !x.range[0]  and x4lim le !x.range[1], xcount )
  if xcount le 0 then begin
     print, ' *** x limits not allowing any data
;;;     return
     !x.range[0]=x4lim[0]
     !x.range[1]=x4lim[ n_elements(x4lim)-1]
  endif

  yinds = where( y4lim ge !y.range[0]  and y4lim le !y.range[1], ycount )
  if ycount le 0 then begin
     print, ' *** y limits not allowing any data
     ;;;return
     !y.range[0]=y4lim[0]
     !y.range[1]=y4lim[ n_elements(y4lim)-1]
  endif
  
  if n_elements( xinds ) gt 1 then begin
     xvalues_opt= xvalues_opt[xinds]
     yvalues_opt= yvalues_opt[yinds]
     x1 = xinds[0]
     x2 = xinds[n_elements(xinds)-1]
     y1 = yinds[0]
     y2 = yinds[n_elements(yinds)-1]

     subIm = image[x1:x2, y1:y2] * 1.0
  endif else begin
     subIm = image
  endelse
  
  ainfo = size( subIm )
  
  maximage = max(subIm, MIN=minimage) 
  if n_elements( bar_range ) eq 2 then begin
     maxbar = bar_range[1]
     minbar = bar_range[0]
  endif else begin
     maxbar = maximage
     minbar = minimage
  endelse
  
  ;; the defaults for the top and bottom colours allow the image colours to be
  ;; independent of the axes and background colours !P.color and !P.background
  ;; (normally 0 and !D.n_colors-1) 
  ;; NB: cant use KEYWORD_SET as 0 is a valid entry !!
  IF N_ELEMENTS(topc) EQ 0 THEN topc =  n_nonfixed-1 ;Brightest color
                                ;(but allow room for plot and axes colours) 
  IF N_ELEMENTS(bottomc) EQ 0 THEN bottomc = 1 ;Dullest color
  
  ;;ensure that the top colour is in the colour table 
  topc = ((topc > 1) <(!D.n_colors -1) ) <(!D.table_size -1) 
  ;;and the bottom colour is less than the top
  bottomc = (bottomc <(topc -1) ) > 0
  ;;help,!d,/st

  sz = size( subIm )                  ;Size of image
  IF sz[0] LT 2 THEN BEGIN
     message, 'Parameter not 2D'
     if (debug) then stop
  ENDIF

  if KEYWORD_SET(debug) then begin
     IF btest(!x.style,1) then  message, /INFO, ' *** !x.style=2 is a problem'
     IF btest(!y.style,1) then  message, /INFO, ' *** !y.style=2 is a problem'
  endif

  ;; NB: cant use KEYWORD_SET as 0 is a valid entry !!
  IF (N_ELEMENTS(bar_ticks) LE 0) THEN bar_ticks = !Z.ticks
  IF (NOT KEYWORD_SET(bar_tickv) ) THEN $
    IF (bar_ticks GT 0) THEN bar_tickv = !Z.tickv(0:bar_ticks) 
  IF (NOT KEYWORD_SET(bar_ticklen) ) THEN bar_ticklen = !Z.ticklen
  IF (NOT KEYWORD_SET(bar_tickname) ) THEN bar_tickname = !Z.tickname

  IF (!X.charsize GT 0) THEN xcharsize = !X.charsize ELSE xcharsize = 1.0
  IF (!P.charsize GT 0) THEN pcharsize = !P.charsize ELSE pcharsize = 1.0
  IF (NOT KEYWORD_SET(tsize) ) THEN tsize = pcharsize
  IF (N_ELEMENTS(raise_ptitle) LE 0 ) THEN raise_ptitle = 1

  pposition = !P.position
  pmulti = !P.multi
  
  IF (KEYWORD_SET(cont) ) THEN nobar = 1
  IF (KEYWORD_SET(nobar) ) THEN put_scale_bar = 0 ELSE put_scale_bar = 1
  
  ;;use plot to set the full window dimensions
  IF (!X.range[0] EQ !X.range[1] ) THEN xtmp = [ 0, sz[1]-1 ] $ ;-0.5 $ 25/10/94
  ELSE xtmp = !X.range
  IF (!Y.range[0] EQ !Y.range[1] ) THEN ytmp = [ 0, sz[2]-1 ] $ ;-0.5 $ 25/10/94
  ELSE ytmp = !Y.range

  IF (KEYWORD_SET(debug) ) THEN BEGIN
     print,'pos 2 !x.range=', !x.range
;;;     print,'pos 2 xrange=', xrange
     print,'pos 2 xtmp=', xtmp
  ENDIF

  plot, xtmp, ytmp, /NODATA, XSTYLE=5, YSTYLE=5,  $
    TITLE=' ', SUBTITLE=' ', XTITLE=' ', YTITLE=' ',  $
	  xtickformat=xtickformat, ytickformat=ytickformat,   $
    NOERASE=KEYWORD_SET(noerase) 
    ;;XRANGE = xtmp, YRANGE=ytmp
                                ;get the full window dimensions
  xwin_sz = !X.window
  ywin_sz = !Y.window
                                ;x-char size in normalised coords is:
  xchar_sz = float(!D.x_ch_size) /float(!D.x_vsize) 
  ychar_sz = float(!D.y_ch_size) /float(!D.y_vsize) 
  
  ;;redefine the y size to allow enough room to put a Plot Title (!P.title) on
  ;;the image when the x-ticks are pointing out ! (this should really be
  ;;handled by the IDL internal routines but isnt YET !)
  yshrink_for_title = 3 * tsize * xchar_sz
  yshrink_for_subtitle = 0
  IF N_ELEMENTS(_extra) GT 0 THEN BEGIN
     dum = WHERE( TAG_NAMES(_extra) EQ 'SUBTITLE', nfound)
     IF nfound GT 0 THEN BEGIN
	subtitle = _extra.SUBTITLE
	IF STRLEN( subtitle ) GT 0 THEN $
	        yshrink_for_subtitle = 2. * tsize * ychar_sz
     ENDIF 
  ENDIF
  
  bar_label_space = 2. * xchar_sz

  IF N_ELEMENTS(bar_label_length) LE 0 THEN BEGIN
     bar_more_space = 0
  ENDIF ELSE BEGIN
     bar_more_space = bar_label_length * xchar_sz
  ENDELSE
  
  !P.position = [ xwin_sz[0],                        $
                  ywin_sz[0] + yshrink_for_subtitle, $
                  xwin_sz[1] - bar_more_space,       $
                  ywin_sz[1] - yshrink_for_title]
  outPos = !P.position 
                                ;size of window for scale bar
                                ;this needs to be outside put_scale_bar loop
                                ;for the recursive call ! 
  IF (N_ELEMENTS(barsz_chars) LE 0) THEN barsz_chars = 2.
  bar_sz = barsz_chars * xchar_sz

  IF (put_scale_bar) THEN BEGIN
    
    ;; NB: cant use KEYWORD_SET as 0 is a valid entry !!
    IF N_ELEMENTS(bar_separation) GT 0 THEN BEGIN
      bar_image_separation = bar_separation * xchar_sz
    ENDIF ELSE BEGIN
      bar_image_separation = 2. * xchar_sz
    ENDELSE

    ;;now redefine the plot window width to allow enough room for a scale bar
    ;;!P.position[2] = !P.position[2] - (bar_sz*4.5) - (-2.*!Y.ticklen > 0.0)
    !P.position[2] = !P.position[2] $ ;the user defined position of RHS
      - abs(bar_sz)  $               ;but allow enough room for the colour bar
      -( -2. * !Y.ticklen > 0.0) $ ;and the ticks on both the bar and image
      - bar_image_separation $  ;and space between bar and image and
      - bar_label_space         ;space for the label on the bar and the title
    
    IF !P.position[2] LT !P.position[0] THEN BEGIN 
      IF (KEYWORD_SET(debug) ) THEN message, /INFO, ' WARNING...'  
      message, /INFO, 'There is NOT enough room for the colour scale bar' 
      message, /INFO, 'to fit into the plotting window !!!' 
      !P.position[2] = !P.position[2] >(!P.position[0] +xchar_sz) 
;;;stop
    ENDIF
    outPos = !P.position 
  ENDIF                         ;put_scale_bar
  
  ;;if bar_tickv set then force the image to have at least this range
  IF (KEYWORD_SET(bar_tickv) ) THEN BEGIN
    maxbar = max([ bar_tickv, maximage, minimage], MIN=minbar) 
    ;;        IF (maxbar GT maximage) THEN a(sz[1]-1, sz[2]-1) = maxbar
    ;;        IF (minbar LT minimage) THEN a(sz[1]-1, 0) = minbar
  ENDIF ELSE BEGIN
    maxbar = maximage
    minbar = minimage
  ENDELSE

  ;;if bar_range set then force the image to have exactly this range
  ;; (over-rides the bar_tickv forcing)
  IF (KEYWORD_SET(bar_range) ) THEN BEGIN
    IF (bar_range[0] NE bar_range[1] ) THEN BEGIN
      maxbar = bar_range[1] & minbar = bar_range[0] 
      ;;            IF (maxbar GT maximage) THEN a(sz[1]-1, sz[2]-1) = maxbar
      ;;            IF (minbar LT minimage) THEN a(sz[1]-1, 0) = minbar
      ;;            ;;clip the image array if necessary
      ;;            a = a > minbar < maxbar
      ;;            ;;the array is now clipped at the bytscl stage

                                ;if (keyword_set(debug)) then begin
      IF (maxbar LT maximage) THEN BEGIN
        IF (KEYWORD_SET(debug) ) THEN BEGIN
          message, /INFO, $
            " WARNING...Maximum image value out of range, Image CLIPPED"
        ENDIF ELSE BEGIN
          message, /INFO, $
            "...Maximum image value out of range, Image CLIPPED"
        ENDELSE
      ENDIF
      IF (minbar GT minimage) THEN BEGIN
        IF (KEYWORD_SET(debug) ) THEN BEGIN
          message, /INFO, $
            " WARNING...Minimum image value out of range, Image CLIPPED"
        ENDIF ELSE BEGIN
          message, /INFO, $
            "...Minimum image value out of range, Image CLIPPED"
        ENDELSE
      ENDIF

    ENDIF                       ;valid bar_range
  ENDIF                         ; bar_range

  IF (KEYWORD_SET(debug) ) THEN BEGIN
     print,'POS 3 !x.range=', !x.range
;;;     print,'POS 3 xrange=', xrange
     print,'pos 3 xtmp=', xtmp
  ENDIF

  ;;use plot to set the scaling with the new window
  comboStruct = ADD_STYLE( _extra, XSTYLE=5, YSTYLE=5 )
  comboStruct.XSTYLE = comboStruct.XSTYLE AND 'FFFD'x	; style=2 not allowed
  comboStruct.YSTYLE = comboStruct.YSTYLE AND 'FFFD'x	; style=2 not allowed
  plot, xtmp, ytmp, /NODATA,  TITLE=' ', SUBTITLE=' ',  $
        XTITLE=' ', YTITLE=' ', /NOERASE,  $
	  xtickformat=xtickformat, ytickformat=ytickformat,   $
	_extra=comboStruct
	
   ;;;help, /str, comboStruct
  ;;XRANGE=xtmp, YRANGE=ytmp

  ;;MOVED to the bytscl code ~20 lines down
  ;; reverse the colour sense of the image if producing greyscale PS output 
  ;;if (!d.name eq 'PS') then a = -a
  
  grey_scale_PS = 0b            ;default is IN COLOUR for ALL devices
  IF (!D.name EQ 'PS') THEN BEGIN ;the output device is PostScript
    IF KEYWORD_SET(debug) THEN message, /INFO, 'PostScript device Detected'
    ;;if the PSETUP suite of routines have been used and the colour variable
    ;;has been set then dont worry about checking for grey-scales
    IF N_ELEMENTS(colour) GT 0 THEN BEGIN
      IF KEYWORD_SET(debug) THEN message, /INFO, '"colour" variable is defined'
      IF KEYWORD_SET(debug) THEN help, colour
      grey_scale_PS = NOT colour
    ENDIF ELSE BEGIN
      ;;decide if the default grey-scale has been loaded
      IF KEYWORD_SET(debug) THEN message, /INFO, 'Getting current colour table'
      tvlct, /GET, r, g, b      ;load in the colour table
      grey_scale_PS = min((r AND g AND b) EQ r) 
    ENDELSE
    IF KEYWORD_SET(debug) THEN help, grey_scale_PS
  ENDIF

  image_window = [ !X.window[0], !Y.window[0], !X.window[1], !Y.window[1] ]

  px = !X.window * !D.x_vsize	;Get size of window in device units
  py = !Y.window * !D.y_vsize
  swx = px[1] -px[0]            ;Size in x in device units
  swx_img = swx
  swy = py[1] -py[0]            ;Size in Y
  swx_offset = 0.0
  if n_elements( xrange ) eq 2 then begin
     if xrange[1] ne xrange[0] and n_elements( xvalues_opt ) gt 0 then begin
        minX = MIN( xvalues_opt, max=maxX )
	swx_img = swx * (maxX-minX )/(xrange[1]-xrange[0])
	if minX gt xrange[0] then swx_offset = float(minX - xrange[0])/  $
	                                       (xrange[1]-xrange[0]) * $
	                                       (px[1] -px[0])
     endif
  endif
  six = float(sz[1] )           ;Image sizes
  siy = float(sz[2] ) 
  aspi = six / siy		;Image aspect ratio
  aspw = swx_img / swy		;Window aspect ratio
  f = aspi / aspw               ;Ratio of aspect ratios

  IF (!D.name NE 'TEK') AND (NOT KEYWORD_SET(cont) ) THEN BEGIN 
    ;;we can use tvscl

;;;print,'before tv'
    IF (KEYWORD_SET(debug) ) THEN message, /INFO, ' .......... Imaging....'

    IF (grey_scale_PS) THEN BEGIN
      IF (KEYWORD_SET(debug) ) THEN BEGIN
        message, /INFO, 'Grey-Scale PostScript Image...'
        message, /INFO, '   ... reversing colour sense of image so dark '
        message, /INFO, '   backgrounds come out light'
      ENDIF
      
      subIm = -subIm                    ;reverse the colour sense of the image if
                                ;producing postscript output so that those
                                ;dark backgrounds come out light !!
      
      ;;rescale into the allowed colour range then shift up to the bottom
      ;;colour index 
      bbb = bytscl(a, MIN=-maxbar, MAX=-minbar, TOP=topc -bottomc) +  $
                   bottomc
      
    ENDIF ELSE BEGIN
      IF (KEYWORD_SET(debug) ) THEN BEGIN
        message, /INFO, 'Colour Image...  '
      ENDIF
      ;;rescale into the allowed colour range then shift up to the bottom
      ;;colour index 
      bbb = bytscl( subIm, MIN=minbar, MAX=maxbar, TOP=topc -bottomc) + bottomc
IF (KEYWORD_SET(debug) ) THEN help,topc,bottomc
      ;;if (!d.name eq 'PS') then bbb(where(bbb eq 1b)) = topc
    ENDELSE
    
    IF (KEYWORD_SET(debug) ) THEN BEGIN
      fstr =  $
        '("[",f8.2,",",f8.2,"] ==> [",f8.2,",",f8.2,"] ==> [",i3,",",i3,"]")'
      message, /INFO, 'byte scaled....' + $
        strcompress( string( minimage, maximage,  $
                             minbar, maxbar, min(bbb), max(bbb), $
                             FORM=fstr) ) 
                                ;help,/dev
    ENDIF
    IF (!D.flags AND 1) NE 0 THEN BEGIN	;Scalable pixels --------------------->
      IF (KEYWORD_SET(debug) ) THEN  $
        message, /INFO, 'Scalable pixels detected'
      
      IF KEYWORD_SET(aspect) THEN BEGIN ;Retain aspect ratio?
        IF f GE 1.0 THEN swy = swy / f  $
        ELSE swx = swx * f      ;Adjust window size
      ENDIF
      
         ; interpolating to the resolution of Postscript files takes 10s of sec
	 ;   and results in huge files, so limit to 2000 pixels in either directions
      
      maxRes = 2000.
      if swX gt maxRes or swY gt MaxRes then shrink = (maxRes/swx) < (maxRes/swy) $
                                        else shrink = 1.0
      fakeSwX = swX*shrink
      fakeSwY = swY*shrink
      newx = findgen(fakeSwX)/fakeSwX*(sz[1]-1)
      newy = findgen(fakeSwY)/fakeSwY*(sz[2]-1)

      resampled_image = INTERPOLATE(subIm, newx, newy, $
			 /GRID, CUBIC=KEYWORD_SET(cubic) ) 
      bbb = bytscl( resampled_image, MIN=minbar, MAX=maxbar, TOP=topc -bottomc) $
              + bottomc
IF (KEYWORD_SET(debug) ) THEN help,topc,bottomc
      tv, bbb, px[0], py[0], XSIZE= SwX, YSIZE= SwY, /DEVICE
      if debug then help,swx,swy,fakeSwX,fakeSwY,newx,newy,subim,bbb
;;;stop
    ENDIF ELSE BEGIN            ;END scalable pixels	<-----------------------
      IF (KEYWORD_SET(debug) ) THEN  $
        message, /INFO, 'NON-Scalable pixels detected'
      IF KEYWORD_SET(window_scale) THEN BEGIN ;Scale window to image?
        tv, bbb, px[0], py[0]   ;Output image
;;;stop
        swx = six		;Set window size from image
        swy = siy
      ENDIF ELSE BEGIN          ;End Scale window to image
        IF (KEYWORD_SET(debug) ) THEN  $
          message, /INFO, 'Scale image to window'
        IF KEYWORD_SET(aspect) THEN BEGIN
          IF f GE 1.0 THEN swy = swy / f ELSE swx = swx * f
        ENDIF                   ;aspect
        
        ;;Have to resample image
        IF KEYWORD_SET(congrid_keyword) THEN BEGIN
;;;          resampled_image = congrid(bbb, $
          resampled_image = congrid(subim, $
                                /Minus_One,  $ ; else will be off by 1/2
                                    swx_img,  $
                                    swy,  $
                                    INTERP=KEYWORD_SET(interp), $
                                    CUBIC=KEYWORD_SET(cubic) )
        ENDIF ELSE BEGIN
	  ; I couldn't get poly_2d to interpolate to the grid points!
;;;          resampled_image = poly_2d(bbb, $	
	  ; BD change 4/00
          IF NOT KEYWORD_SET(interp) THEN BEGIN
             resampled_image = poly_2d(subIm, $
                                       [[ 0, 0], [ six /swx_img, 0]],  $
                                       [[ 0, siy /swy], [ 0, 0]], $
                                       KEYWORD_SET(interp),  $
                                       (swx_img>2),  $
                                       (swy>2),  $
                                       CUBIC=KEYWORD_SET(cubic) )
          ENDIF ELSE BEGIN
	     newx = findgen(swx)/swx*(sz[1]-1)
	     newy = findgen(swy)/swy*(sz[2]-1)
             resampled_image = INTERPOLATE(subIm, newx, newy, $
                                	   /GRID, CUBIC=KEYWORD_SET(cubic) )
          ENDELSE
        ENDELSE
	
	; BD change 4/00
	subIm = congrid(resampled_image,sz[1],sz[2],INTERP=KEYWORD_SET(interp))	; so contours match colors
        reim_sc = bytscl(resampled_image, MIN=minbar, MAX=maxbar, $
	                 TOP=topc -bottomc) + bottomc
        IF (KEYWORD_SET(debug) ) THEN help,topc,bottomc

        tv, reim_sc, px[0]+1+swx_offset, py[0]+1 ; BD change 12/6/99
        
      ENDELSE			;End Scale image to window
    ENDELSE                     ;scalable pixels
  ENDIF

  ;;use contour to give the axes all the time
  ppmulti = !P.multi
  !P.multi = pmulti
  ; BD change 4/00 because bug if 10x10, eg.,
  xstyle = 1
  IF ((!X.style MOD 2) EQ 0) THEN xstyle = 1 +!X.style ELSE xstyle = !X.style
  IF (!Y.style EQ 2) THEN ystyle = 1
  ; BD change 4/00 because bug if 10x10, eg.,
  ystyle=1
  IF ((!Y.style MOD 2) EQ 0) THEN ystyle = 1 +!Y.style ELSE ystyle = !Y.style
  c_sz = sz

     if (debug) then print,'pos 4 xtmp=', xtmp

  IF KEYWORD_SET(exact) THEN BEGIN
    c_x = ((indgen(c_sz[1] ) +0.5) *(xtmp[1] -xtmp[0] ) /(c_sz[1] ) )  $
      +xtmp[0] 
    c_y = ((indgen(c_sz[2] ) +0.5) *(ytmp[1] -ytmp[0] ) /(c_sz[2] ) )  $
      +ytmp[0] 
  ENDIF ELSE BEGIN
    c_x = (indgen(c_sz[1] ) *(xtmp[1] -xtmp[0] ) /(c_sz[1] -1 ) ) +xtmp[0] 
    c_y = (indgen(c_sz[2] ) *(ytmp[1] -ytmp[0] ) /(c_sz[2] -1 ) ) +ytmp[0] 
  ENDELSE
  
  IF Use_Args THEN BEGIN
     c_x = xvalues_opt
     c_y = yvalues_opt
  ENDIF
  
  comboStruct = ADD_STYLE( _extra, XSTYLE=xstyle, YSTYLE=ystyle )
  comboStruct.XSTYLE = comboStruct.XSTYLE AND 'FFFD'x	; style=2 not allowed
  comboStruct.YSTYLE = comboStruct.YSTYLE AND 'FFFD'x	; style=2 not allowed
  IF (KEYWORD_SET(debug) ) THEN BEGIN
     print,'!x.crange=', !x.crange
     print,' '
  ENDIF
  
  ;;;help, a, c_x, c_y
  
  contour, subIm, c_x, c_y, /NOERASE, /NODATA,  $
       XRANGE=!x.crange, YRANGE=!y.crange, $       ; not sure why needed
       POS=[ px[0]/!d.x_size, py[0]/!d.y_size, $
            (px[0]+swx)/!d.x_size, (py[0]+swy)/!d.y_size], $
	  xtickformat=xtickformat, ytickformat=ytickformat,   $
       TITLE=' ', SUBTITLE=' ', _EXTRA=comboStruct
;;;stop
;;;print, 'c_x at contour 1='
;;;minmax, c_x

  ;;find the position for the title and put it there
  xpos = 0.5 * total(!X.window) 
  ;;if (!x.charsize gt 0) then xcharsize=!x.charsize else xcharsize=1.0
  ;;if (!p.charsize gt 0) then pcharsize=!p.charsize else pcharsize=1.0
  ;;if (not keyword_set(tsize)) then tsize = pcharsize
  ch_dev_to_norm = pcharsize * !D.y_ch_size /!D.y_vsize 
  ;;ypos = !y.window[1] + (-!x.ticklen > 0.0) +2.7*xcharsize*ch_dev_to_norm
  ypos = !Y.window[1] +( -!X.ticklen > 0.0)  $
    +raise_ptitle * xcharsize * ch_dev_to_norm
  xyouts, xpos, ypos, !P.title, ALIGNMENT=0.5, /NORM, SIZE=tsize * 1.2

  ;;find the position for the subtitle and put it there
  ypos = !Y.window[0] -( -!X.ticklen > 0.0)  $
    -(2.7 * xcharsize +tsize) * ch_dev_to_norm
  xyouts, xpos, ypos, !P.subtitle, ALIGNMENT=0.5, /NORM, SIZE=tsize

  cloop = 0
  IF (NOT KEYWORD_SET(nocont) ) THEN BEGIN

    ;;use contour to give data contours on top of the image
    IF (KEYWORD_SET(debug) ) THEN message, /INFO, ' .......... Contouring....'

    ;;need to reset the ranges because the axis labels are independent 
    ;;of the way contour works (plots as index numbers)
    save_p = !P
    save_x = !X
    save_y = !Y	
    save_z = !Z	
    resetplt, /ALL
    clearplt, /ALL
    !X.style = 5
    !Y.style = 5
    !P.multi = pmulti
    IF (!P.thick LE 0) THEN !P.thick = 1

                                ;plot contours
    IF (KEYWORD_SET(levels) ) THEN csz = N_ELEMENTS(levels) $
    ELSE IF (KEYWORD_SET(nlevels) ) THEN csz = nlevels +1 ELSE csz = 6
    IF (NOT KEYWORD_SET(max_value) ) THEN max_value = maximage +1

    keyc_colors = KEYWORD_SET(c_colors) 
    keyc_linestyle = KEYWORD_SET(c_linestyle) 
    keyc_thick = KEYWORD_SET(c_thick) 

    cont_image = image * 1.0
    IF (KEYWORD_SET(badpts) ) THEN BEGIN
      IF (badpts[0] GE 0) THEN cont_image(badpts) = 10 * maximage
    ENDIF
    
        if (debug) then  print,'pos 5 xtmp=', xtmp

    c_sz = size(cont_image) 
    IF KEYWORD_SET(exact) THEN BEGIN
      c_x = ((indgen(c_sz[1] ) +0.5) *(xtmp[1] -xtmp[0] ) /(c_sz[1] ) )  $
        +xtmp[0] 
      c_y = ((indgen(c_sz[2] ) +0.5) *(ytmp[1] -ytmp[0] ) /(c_sz[2] ) )  $
        +ytmp[0] 
    ENDIF ELSE BEGIN
      c_x = (indgen(c_sz[1] ) *(xtmp[1] -xtmp[0] ) /(c_sz[1] -1 ) ) +xtmp[0] 
      c_y = (indgen(c_sz[2] ) *(ytmp[1] -ytmp[0] ) /(c_sz[2] -1 ) ) +ytmp[0] 
    ENDELSE
      
  
    IF Use_Args THEN BEGIN
       c_x = xvalues_opt
       c_y = yvalues_opt
    ENDIF

    IF (KEYWORD_SET(levels) ) THEN BEGIN
      cloop = 1
      IF (NOT keyc_colors) THEN c_colors =  $
        (fix((((levels -minbar) /float(maxbar -minbar) ) * 0.5 + 0.75)  $
             *(topc -bottomc) ) MOD(topc -bottomc +1) ) + bottomc
      IF (NOT keyc_linestyle) THEN c_linestyle = levels * 0
      IF (NOT keyc_thick) THEN c_thick = levels * 0 +!P.thick
      zero = where(levels EQ 0, count) 
      zero = zero[0] 
      IF (count * zero GT 0) THEN BEGIN
        ;; do special things to highlight the zero contour
        IF (NOT keyc_colors) AND (KEYWORD_SET(cont) ) THEN $
          c_colors(zero) = topc
        IF (NOT keyc_thick) THEN BEGIN
          IF (KEYWORD_SET(cont) ) THEN  $
            c_thick(zero) = 2.5 * !P.thick  $
          ELSE  $
            c_thick(zero) = 2.0 * !P.thick
        ENDIF
        IF (NOT keyc_linestyle) AND  $
          (KEYWORD_SET(cont) ) THEN c_linestyle(0:zero -1) = 1 
      ENDIF
      IF (NOT keyc_colors) AND (grey_scale_PS) AND (KEYWORD_SET(cont) ) THEN  $
        c_colors = topc
                                ;invert line colors if b/w postscript output
      IF (grey_scale_PS) THEN usec_cols = topc - c_colors + bottomc $
      ELSE usec_cols = c_colors
      comboStruct = ADD_STYLE( _extra, XSTYLE=xstyle, YSTYLE=ystyle )
      comboStruct.XSTYLE = comboStruct.XSTYLE AND 'FFFD'x	; style=2 not allowed
      comboStruct.YSTYLE = comboStruct.YSTYLE AND 'FFFD'x	; style=2 not allowed

      contour, cont_image, c_x, c_y, /overplot,  $
        XRANGE=!x.crange, YRANGE=!y.crange, $   ; BD add 4/00
        POS=[ px[0], py[0], px[0] +swx, py[0] +swy], $
        MAX_VALUE=max_value, C_COLORS=usec_cols, $
        C_LINESTYLE=c_linestyle, C_THICK=c_thick, $
	  xtickformat=xtickformat, ytickformat=ytickformat,   $
        LEVELS=levels, _EXTRA=comboStruct
;;;help,cont_image, c_x, c_y
;;;print, 'c_x:'
;;;minmax,c_x
;;;stop


    ENDIF ELSE BEGIN
      cloop = 2
      IF (NOT keyc_colors) THEN c_colors =  $
        [ replicate(topc, (csz /2 > 1) ),  $
          replicate(bottomc +(topc -bottomc) /4, (csz /2 +1 > 1) ) ]
      IF (NOT keyc_linestyle) THEN c_linestyle = replicate(0, csz) 
      IF (NOT keyc_thick) THEN c_thick = [ !P.thick]
      IF (NOT KEYWORD_SET(nlevels) ) THEN nlevels = csz -1
      
                                ;invert line colors if b/w postscript output
      IF (grey_scale_PS) THEN usec_cols = topc - c_colors + bottomc $
      ELSE usec_cols = c_colors

;;;stop
      comboStruct = ADD_STYLE( _extra, XSTYLE=xstyle, YSTYLE=ystyle )
      comboStruct.XSTYLE = comboStruct.XSTYLE AND 'FFFD'x	; style=2 not allowed
      comboStruct.YSTYLE = comboStruct.YSTYLE AND 'FFFD'x	; style=2 not allowed
      contour, cont_image, c_x, c_y, /overplot, $
        XRANGE=!x.crange, YRANGE=!y.crange, $   ; BD add 4/00
        MAX_VALUE=max_value, C_COLORS=usec_cols, $
        C_LINESTYLE=c_linestyle, C_THICK=c_thick, $
	  xtickformat=xtickformat, ytickformat=ytickformat,   $
        NLEVELS=nlevels, _EXTRA=comboStruct
;;;help,cont_image, c_x, c_y
;;;print, ' from way down, c_x:'
;;;minmax,c_x

;;;stop


    ENDELSE

if n_elements( currentTime ) gt 0 then $
   oplot, [currentTime,currentTime], !y.crange, line=2, thick=2

    cont_image = 0              ;release some space NOW

    !P = save_p
    !X = save_x
    !Y = save_y
    !Z = save_z
  ENDIF

if n_elements( currentTime ) gt 0 then $
   oplot, [currentTime,currentTime], !y.crange, line=2, thick=2

  IF put_scale_bar AND (bar_sz GT 0) THEN BEGIN
    ;;add a scale bar by calling this routine AGAIN
    
    IF (KEYWORD_SET(debug) ) THEN  $
      message, /INFO, ' .......... adding Scale Bar....'

    ;;do it !!!
    ;;now redefine the plot position for a scale bar
    ;;Note that the y pos must be re-enlarged for the shrinkage that took place
    ;;for the Plot title, as the recursive call will do the shrinkage AGAIN !!
    
    !P.position = [ !X.window[1] +  $
                    bar_image_separation + ( -!Y.ticklen > 0.0), $
                    !Y.window[0], $
                    !X.window[1] +  $
                    abs(bar_sz) + bar_label_space +( -!Y.ticklen > 0.0), $
                    !Y.window[1] + yshrink_for_title]
                                ;use image_cont again to do scale bar
    save_p = !P
    save_x = !X
    save_y = !Y
    clearplt, /X, /Y, /P
    !P.multi = pmulti
    !P.thick = save_p.thick
    !P.ticklen = save_p.ticklen
    !X.minor = -1 & !X.ticklen = 0.000001
                                ;if bar_tickv set then force the image to
                                ;have at least this range
                                ;maxscl = max([maxbar,maximage])
                                ;minscl = min([minbar,minimage])
    ;;bardat = (findgen(!D.n_colors) /(!D.n_colors -1) *(maxbar -minbar) ) $
    ;;  + minbar
    bar_dim = !D.n_colors < !D.table_size
    bardat = (findgen(bar_dim) /(bar_dim -1) *(maxbar -minbar) ) + minbar
    bardat = transpose([[ bardat], [ bardat]]) 
    !X.style = 1
    !Y.style = 1
    
    !P = save_p
    CASE (cloop) OF
      
      0: th_image_cont_4bar, bardat, /NOBAR, /NOERASE,  $
        NOCONT=KEYWORD_SET(nocont), TOPC=topc, BOTTOMC=bottomc
      
      1: th_image_cont_4bar, bardat, /NOBAR, /NOERASE,  $
        NOCONT=KEYWORD_SET(nocont), TOPC=topc, BOTTOMC=bottomc, $
        MAX_VALUE=max_value, C_COLORS=c_colors, $
        C_LINESTYLE=c_linestyle, C_THICK=c_thick, LEVELS=levels
      
      2: th_image_cont_4bar, bardat, /NOBAR, /NOERASE,  $
        NOCONT=KEYWORD_SET(nocont), TOPC=topc, BOTTOMC=bottomc, $
        MAX_VALUE=max_value, C_COLORS=c_colors, $
        C_LINESTYLE=c_linestyle, C_THICK=c_thick, NLEVELS=nlevels
      
    ENDCASE
    clearplt,/p
    resetplt, /Y
    !P.charsize = save_p.charsize
    !Y.charsize = save_y.charsize
    !Y.ticklen = save_y.ticklen
    IF (N_ELEMENTS(bar_ticks) GT 0) THEN !Y.ticks = bar_ticks ELSE $
      IF (NOT KEYWORD_SET(bar_range) ) AND ((maxbar -minbar) LE 10.0) $
      THEN !Y.ticks = 2
    IF (KEYWORD_SET(bar_tickv) ) THEN BEGIN
      !Y.tickv = bar_tickv
      !Y.ticks = N_ELEMENTS(bar_tickv) -1
    ENDIF
    !Y.range = [ minbar, maxbar]
    IF (KEYWORD_SET(bar_tickname) ) THEN !Y.tickname = bar_tickname
    IF (KEYWORD_SET(bar_ticklen) ) THEN !Y.ticklen = bar_ticklen
    !Y.ticklen = !Y.ticklen *  $
      (image_window[2] -image_window[0] )  $
      /(save_p.position[2] -save_p.position[0] ) 
;;;    black = mk_color('black',/search)
    axis, /YAX, YTITLE=!Z.title, YSTYLE=1, /NOERASE, ytickformat='betterticklabels', $
          color=!p.color
;;;          color=black
    
    !P = save_p
    !X = save_x
    !Y = save_y
  ENDIF

   !P = orig_p
   !X = orig_x
   !Y = orig_y       
   !Z = orig_z

  !P.position = pposition
  !P.multi = ppmulti
  
  RETURN
END