;+
; NAME:
;   plotsub
;
; PURPOSE:
;   plot a subarray with pixel axes shown
;
; CATEGORY:
;   Plotting
;
; CALLING SEQUENCE:
;   plotsub, d, xm, ym, wx, wy, bin, magnify=magnify, $
;             windx=windx, zmin=zmin, zmax=zmax, fullscale=fullscale, $
;             dsub=dsub, title=title, autosub=autdosub
;
;
; INPUTS:
;       D:      The array to be plotted
;       XM:     Center of subarray to be plotted
;       YM:     Center of subarray to be plotted
;       WX:     Full width of subarray to be plotted
;       WY:     Full width of subarray to be plotted
;       BIN:    Magnification or binning factor (left here for compat.)
;
; OPTIONAL KEYWORD PARAMETERS:
;  magnify - 1 if magnifying
;  windx - window index (otherwise, current window)
;  zmin, zmax - min mad max to plot
;  fullscale - set zmin and zmax to min and max of image
;  title - title for window
;  bin - Magnification or binning factor
;  autosub - region to use for scaling
;
; OPTIONAL OUTPUT PARAMETERS:
;  dsub - data subarray
;  pos - plot position
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None.
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;   2006 Mar 20 Leslie Young SwRI
;   2012 May 2 LAY retain=2 for window
;-

pro plotsub, d, xm, ym, wx, wy, binp, bin=bin, magnify=magnify, $
             windx=windx, zmin=zmin, zmax=zmax, fullscale=fullscale, $
             invert=invert, $
             dsub=dsub, pos=pos, $
             title=title, notick=notick, autosub=autosub, dx0=dx0

    device, decomp=1

    ; is this a 2-D gray-scale or a 3-D true color?
    ndim = size(d,/n_dim)
    dim = size(d, /dim)
    if ndim eq 3 then begin
       if dim[0] eq 3 then begin
          true = 1
          nx_full = dim[1]
          ny_full = dim[2]
       endif else begin
          if dim[2] eq 3 then begin
             true = 3
             nx_full = dim[0]
             ny_full = dim[1]
          endif else begin
             if dim[1] eq 3 then begin
                true = 2
                nx_full = dim[0]
                ny_full = dim[2]
             endif else begin
                print, 'plotsub: d is not 2-D or a true-color 3-D array'
                return
             endelse
          endelse
       endelse
    endif else begin
       true = 0
       nx_full = dim[0]
       ny_full = dim[1]
    endelse

    case  n_params() of
        1: begin
            wx = nx_full
            xm = (wx - 1)/2.
            wy = ny_full
            ym = (wy - 1)/2.
            if not keyword_set(bin) then begin
                bin = max(ceil(dim/(get_screen_size()-[0,200.])))
            end
        end
        5: begin
            if not keyword_set(bin) then begin
                bin = max(ceil([wx,wy]/(get_screen_size()-[0,200.])))
            end
        end
        6: begin
            bin = binp
        end
    end

    if isarray(xm) then xm = xm[0]    
    if isarray(ym) then ym = ym[0]

    ; nominal edges
    if keyword_set(magnify) then begin
        wxbin = floor(wx*bin)   ; number of binned pixels
        wybin = floor(wy*bin) 
        hx = (wxbin-1)/2./bin   ; half width
        hy = (wybin-1)/2./bin   ; half width
    endif else begin
        wxbin = floor(wx/bin)   ; number of binned pixels
        wybin = floor(wy/bin) 
        hx = bin * (wxbin-1)/2. ; half width
        hy = bin * (wybin-1)/2. ; half width
    endelse
    x0 = round(xm - hx)
    y0 = round(ym - hy)
    if keyword_set(magnify) then begin
        x1 = x0 + wxbin/bin - 1
        y1 = y0 + wybin/bin - 1
    endif else begin
        x1 = x0 + bin*wxbin - 1
        y1 = y0 + bin*wybin - 1
    endelse

    ; make sure we're in bounds
    if keyword_set(magnify) then begin
        if x0 lt 0       then x0 = 0
        if x1 ge nx_full then x1 = nx_full-1
        if y0 lt 0       then y0 = 0
        if y1 ge ny_full then y1 = ny_full-1
    endif else begin
        while x0 lt 0       do x0 = x0 + bin
        while x1 ge nx_full do x1 = x1 - bin
        while y0 lt 0       do y0 = y0 + bin
        while y1 ge ny_full do y1 = y1 - bin
    endelse
    nx = 1 + x1 - x0
    ny = 1 + y1 - y0
    if keyword_set(magnify) then begin
        nxbin = nx*bin
        nybin = ny*bin
    endif else begin
        nxbin = nx/bin
        nybin = ny/bin
    endelse
    case true of
       0: dsub = d[x0:x1, y0:y1]
       1: dsub = d[*,x0:x1, y0:y1]
       2: dsub = d[x0:x1, *,y0:y1]
       3: dsub = d[x0:x1, y0:y1,*]
    endcase

    if true ne 0 then begin
       zmin = 0
       zmax = 255
    endif
    if n_elements(zmin) eq 0 or n_elements(zmax) eq 0 then begin
        if keyword_set(fullscale) then begin
            minmax = minmax(dsub)
            if keyword_set(invert) then begin
                zmin = minmax[1]
                zmax = minmax[0]
            endif else begin
                zmin = minmax[0]
                zmax = minmax[1]
            endelse
        endif else begin
            if not keyword_set(autosub) then autosub = dsub
            autoscale,autosub,zzmin,zzmax
            if keyword_set(invert) then begin
                if 0 eq n_elements(zmin) then zmin = zzmax
                if 0 eq n_elements(zmax) then zmax = zzmin
            endif else begin
                if 0 eq n_elements(zmin) then zmin = zzmin
                if 0 eq n_elements(zmax) then zmax = zzmax
            endelse
        endelse
    endif 

    zzmin = zmin < zmax
    zzmax = zmin > zmax
    if keyword_set(invert) then begin
        zmin = zzmax
        zmax = zzmin
    endif else begin
        zmin = zzmin
        zmax = zzmax
    endelse
    
    if keyword_set(magnify) then begin
       case true of
          0: dbin = byte(rebin(linscl(dsub, zmin,zmax,0,254.5, /clip), $
                          nxbin, nybin,/sample) )
          1: dbin = byte(rebin(linscl(dsub, zmin,zmax,0,254.5, /clip), $
                          3, nxbin, nybin,/sample) )
          2: dbin = byte(rebin(linscl(dsub, zmin,zmax,0,254.5, /clip), $
                          nxbin, 3, nybin,/sample) )
          3: dbin = byte(rebin(linscl(dsub, zmin,zmax,0,254.5, /clip), $
                          nxbin, nybin,3, /sample) )
       endcase
    endif else begin
       case true of 
          0: dbin = byte(linscl(rebin(dsub, nxbin, nybin),$
                           zmin,zmax,0,254.5, /clip))
          1: dbin = byte(linscl(rebin(dsub, 3, nxbin, nybin),$
                          zmin,zmax,0,254.5, /clip))
          2: dbin = byte(linscl(rebin(dsub, nxbin, 3, nybin),$
                           zmin,zmax,0,254.5, /clip))
          3: dbin = byte(linscl(rebin(dsub, nxbin, nybin, 3),$
                           zmin,zmax,0,254.5, /clip))
       endcase
    endelse

    ; calculate the margins
    ticklenpix = 10
    x_ch = (!p.charsize>1) * !d.x_ch_size 
    y_ch = (!p.charsize>1) * !d.y_ch_size 
    marg = 10
    if n_elements(dx0) eq 0 then dx0 = marg + ticklenpix + x_ch * ceil(alog10(y1)) + y_ch 
    dy0 = marg + ticklenpix + y_ch * 2
    dx1 = marg + ticklenpix + x_ch * ceil(alog10(x1))/2.
    dy1 = marg + ticklenpix

    pos = round([dx0, dy0, dx0+nxbin-1, dy0+nybin-1])
    xs = pos[2] + dx1
    ys = pos[3] + dy1

    xticklen = -ticklenpix/float(ys)
    yticklen = -ticklenpix/float(xs)

    if n_elements(windx) eq 0 then windx = !d.window
    if windx eq -1 then windx = 0
    if not keyword_set(title) then begin
        if !d.x_size lt xs or !d.y_size lt ys or windx ne !d.window then begin
            window, windx, xs=xs, ys=ys, title=title, retain=2
        endif 
    endif else begin
        window, windx, xs=xs, ys=ys, title=title, retain=2
    endelse

    !p.color = 0
    !p.background = !d.n_colors - 1
    erase
    if keyword_set(magnify) then begin
        xr=[x0-0.5-0.5/bin,x1+0.5+0.5/bin]
        yr=[y0-0.5-0.5/bin,y1+0.5+0.5/bin]
    endif else begin
        xr=[x0-0.5-bin/2.,x1+0.5+bin/2.]
        yr=[y0-0.5-bin/2.,y1+0.5+bin/2.]
    endelse
    if not keyword_set(notick) then begin 
       pos = pos+[-1,-1,1,1]
        plot, [0],[0], xr=xr, yr=yr, $
          pos = pos, /device, /nodata, $
          xticklen=xticklen, yticklen=yticklen, /xs, /ys
    endif else begin
       pos = pos+[-1,-1,1,1]
        plot, [0],[0], xr=xr, yr=yr, $
          pos = pos, /device, /nodata, $
          xticklen=0, yticklen=0, /xs, /ys, $
          xticks=1, yticks=1, xtickname=[' ',' '],ytickname=[' ',' ']
    endelse
    ; at this point, 
    ; pos[0] = !p.clip[0] + 2
    ; pos[1] = !p.clip[1] + 1
    ; pos[2] = !p.clip[2]  
    ; pos[3] = !p.clip[3] - 1

    tv, dbin, pos[0], pos[1], true=true
end