;+
; 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
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None.
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;   2006 Mar 20 Leslie Young SwRI
;-

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

    device, decomp=1

    dim = size(d, /dim)

    case  n_params() of
        1: begin
            wx = dim[0]
            xm = (wx - 1)/2.
            wy = dim[1]
            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 dim[0] then x1 = dim[0]-1
        if y0 lt 0      then y0 = 0
        if y1 ge dim[1] then y1 = dim[1]-1
    endif else begin
        while x0 lt 0      do x0 = x0 + bin
        while x1 ge dim[0] do x1 = x1 - bin
        while y0 lt 0      do y0 = y0 + bin
        while y1 ge dim[1] 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
    dsub = d[x0:x1, y0:y1]

    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
        dbin = byte(rebin(linscl(dsub, zmin,zmax,0,254.5, /clip), $
                          nxbin, nybin,/sample) )
    endif else begin
        dbin = byte(linscl(rebin(dsub, nxbin, nybin),$
                           zmin,zmax,0,254.5, /clip))
    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
    dx0 = marg + ticklenpix + x_ch * ceil(alog10(y1))
    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 
        endif 
    endif else begin
        window, windx, xs=xs, ys=ys, title=title 
    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 
        plot, [0],[0], xr=xr, yr=yr, $
          pos = pos+[-1,-1,1,1], /device, /nodata, $
          xticklen=xticklen, yticklen=yticklen, /xs, /ys
    endif else begin
        plot, [0],[0], xr=xr, yr=yr, $
          pos = pos+[-1,-1,1,1], /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]
end