function Imcut, image, imsize, xc, yc, xc1, yc1, $
                help=help, fixsize=fixsize, $
                peak = peak, $
                silent=silent, pad=pad, mask=mask,$
				center=center,newcenter=newcenter,$
				method2=method2, $
				x0=x0,y0=y0

;+
; PURPOSE:
; Function to cut out a square subsection out of an image,
; with the user specifying the center of the subsection.
;
; NOTES: 
; Will return a non-square image if the specified center is
; too close to the edges.
;
;
; function Imcut, image, imsize, xc, yc, xc1, yc1, $
;                help=help, fixsize=fixsize, $
;                peak = peak, $
;                silent=silent 
;
;
; INPUTS
;   image	original image
;
; OPTIONAL INPUTS
;   imsize	size of the subsection (256 default)
;               can be a scalar or a 2-element vector (x,y)
;   xc,yc	center location of the subsection 
;		(default is center)
; 
; OPTIONAL OUTPUTS
;  xc1,yc1      if too close to edge, returns the coords
;               of the original point (xc,yc) in the new subarray
;  x0,y0		coordinates of the corner of the new subarray
;
; KEYWORD PARAMETERS
;   /fixsize    will adjust center so return an image of desired size
;               (by default it trims the output image size keeping
;                the desired center fixed) 
;
;   /pad		instead of just trimming an image near the edge, pad it
;   			with zeros to get it to the correct size. Return a mask
;   			of good (=1) and bad (=0) pixels in the 'mask' variable.
;   /peak       use peak pixel value for central coordinate
;               then no need to pass (xc,yc)
;
; RETURNS
;	the image subsection
;
; HISTORY
; Written by MCL(UCB): 10/01/95
; 02/10/96 (MCL): changed order of parameters (imsize,xc,yc)
; 04/23/96 (MCL): allowed imcut to be 2d 
; 05/04/96 (MCL): added /fixsize
; 03/04/98 (MCL): returns (xc1,yc1), useful for image labelling
;   -> working?
; 10/20/00 (MCL): gives warning if uses default values for imsize,xc,yc
; 06/11/01 (MCL): added /peak
; 2004-10-21  MDP: Added /pad and mask keywords. 
; 2006-04-20 MDP: Made xc1 and yc1 actually work properly.
;
; -> need to check if handling odd numbered sizes properly
;-

;on_error, 2

if keyword_set(help) or n_params() eq 0  $
  or n_elements(image) eq 0 then begin
	message,'function imcut(image,{imsize},{xc},{yc},[help],[silent])'
endif

sz = size(image)
if n_elements(imsize) eq 0 then begin
    nx =  256
    ny =  256
    message, 'imsize not set, using default of 256', /info
endif else if n_elements(imsize) eq 1 then begin
    nx =  imsize
    ny =  imsize
endif else if n_elements(imsize) eq 2 then begin
    nx = imsize(0)
    ny = imsize(1)
endif else begin
    message, 'imsize must contain 1 or 2 elements!', /info
    stop
endelse

if keyword_set(center) then begin
	xc=center[0]
	yc=center[1]
endif

; center on peak pixel if desired
if keyword_set(peak) then  $
  whereismax, image, xc, yc, /silent


if n_elements(xc) eq 0 then begin
    message, 'x-position not set, choosing midpoint', /info
    xc = sz(1)/2
endif 
if n_elements(yc) eq 0 then begin
    message, 'y-position not set, choosing midpoint', /info
    yc = sz(2)/2
endif

print,"Using center:"+printcoo(xc,yc)

dx = fix(nx/2)
dy = fix(ny/2)
x0 = round(xc) - dx > 0
x1 = round(xc) + dx < (sz(1)-1) 
y0 = round(yc) - dy > 0
y1 = round(yc) + dy < (sz(2)-1)
if not(odd(nx)) then x1 = x1-1
if not(odd(ny)) then y1 =  y1 -1

; revised 2006-04-28 by Marshall
; to better handle getting out odd arrays 
; centered on a half pixel
if keyword_set(method2) then begin
dx = (nx/2.)
dy = (ny/2.)
x0 = round(xc - dx) > 0
x1 = round(xc + dx-1) < (sz(1)-1) 
y0 = round(yc - dy) > 0
y1 = round(yc + dy-1) < (sz(2)-1)
	;stop
endif


if (x1-x0+1) lt nx or (y1-y0+1) lt ny then begin
;	print,imsize,dd
;	print,x1,x0,x1-x0+1
;	print,y1,y0,y1-y0+1
    if not(keyword_set(silent)) then begin
        message, 'x-coords = '+printcoo(x0, x1)+$
          'y-coords = '+printcoo(y0, y1), /info
        message, '* too close to edge: '+strc(x1-x0+1)+$

          ','+strc(y1-y0+1)+' *', /info
    endif
    if keyword_set(fixsize) then begin
        if (round(xc)-dx) lt 0 then begin
            x0 = 0
            x1 =  2*dx-1
        endif
        if (round(xc)+dx) gt (sz(1)-1) then begin
            x0 =  sz(1)-2*dx
            x1 = sz(1)-1
        endif
        if (round(yc)-dy) lt 0 then begin
            y0 = 0
            y1 =  2*dx-1
        endif
        if (round(yc)+dy) gt (sz(2)-1) then begin
            y0 =  sz(2)-2*dy
            y1 = sz(2)-1
        endif
        if not(keyword_set(silent)) then $
          print, '* fixing size: x=', printcoo(x0, x1), ', y=', $
          printcoo(y0, y1) 
    endif 
endif

if not(keyword_set(silent)) then $
    print, 'cutting = ', printcoo(x0, y0), ' to ', printcoo(x1, y1)

if keyword_set(pad) then begin
	padded = fltarr(nx,ny)
	mask = intarr(nx,ny)
	if (round(xc) - dx) lt 0 then  xoffset = (-1)*(round(xc) - dx) else xoffset=0
	if (round(yc) - dy) lt 0 then  yoffset = (-1)*(round(yc) - dy) else yoffset=0

	if x1 lt 0 or y1 lt 0 then return,padded ; nothing good at all!
	sz = size(image(x0:x1,y0:y1))
	padded[xoffset:xoffset+sz[1]-1,yoffset:yoffset+sz[2]-1] = image(x0:x1,y0:y1)
	mask[xoffset:xoffset+sz[1]-1,yoffset:yoffset+sz[2]-1]   = 1
	return,padded
	

endif

xc1 = xc - x0
yc1 = yc - y0
newcenter=[xc1,yc1]

sz = size(image)
if sz[0] gt 2 then return,image[x0:x1,y0:y1,*] $
else return,image[x0:x1,y0:y1]

end