;+
; NAME: cube_rebin 
; PURPOSE:
; 	Rebin a cube in the X and Y dimensions, while leaving the Z dimension
; 	unchanged.
;
; INPUTS:
; KEYWORDS:
; OUTPUTS:
;
; HISTORY:
; 	Began 2007-07-16 17:40:19 by Marshall Perrin 
;-


FUNCTION cube_rebin, cubestruct, rebin=rebin_factor, zrange=zrange

	waveaxis = cube_getwaveaxis(cubestruct)
	sz = size(cubestruct.cube)
	szsmall = sz/rebin_factor & szsmall[waveaxis+1] = sz[waveaxis+1] ; don't resamp wavelen

	ratio = float(szsmall[1:3])/sz[1:3]

	; rebin can only handle INTEGER amounts of resizing, so cut out
	; a subarray which CAN be resampled appropriately
	;  NOTE: this ALWAYS crops the top right part off preferentially.

	subscripts = fltarr(3,2)
	subscripts[*,1] = floor(1./ratio) * szsmall[1:3]-1
	ratio = 1./floor(1./ratio) ; integer version
	if keyword_set(zrange) then begin
		subscripts[waveaxis,*] = zrange 
		szsmall[waveaxis+1] = zrange[1]-zrange[0]+1
	endif

	;---- perform some rebinnning

	newcube = rebin(cubestruct.cube[subscripts[0,0]:subscripts[0,1],subscripts[1,0]:subscripts[1,1],$
		subscripts[2,0]:subscripts[2,1]], $
		szsmall[1], szsmall[2], szsmall[3])

;;	if ~(keyword_set(zrange)) then begin
;;		; Rebin the entire cube
;;
;;		newcube = rebin(cubestruct.cube, szsmall[1], szsmall[2], szsmall[3])
;;
;;
;;	endif else begin 
;;		; rebin only a selected wavelength range.
;;		nwave = zrange[1]-zrange[0]+1 ; inclusive
;;		szsmall[waveaxis] = nwave
;;		case waveaxis of 
;;		1: begin
;;			newcube = rebin(cubestruct.cube[zrange[0]:zrange[1],*,*], szsmall[1], szsmall[2], szsmall[3])
;;		end
;;		2: begin
;;			newcube = rebin(cubestruct.cube[*,zrange[0]:zrange[1],*], szsmall[1], szsmall[2], szsmall[3])
;;		end
;;		3: begin
;;			newcube = rebin(cubestruct.cube[*,*,zrange[0]:zrange[1]], szsmall[1], szsmall[2], szsmall[3])
;;		end
;;		endcase
;;	endelse

	;---- update the FITS header astrometry
	;     Code based on hrebin from Goddard
	extast3, cubestruct.header, astr, numparams

	newastr = astr
	; TODO fix this to be more accurate with what happens up above!
	for ia=0L,2 do begin
		if keyword_set(exact) and (not keyword_set(SAMPLE)) and (ratio[ia] GT 1) then $
     	 newastr.crpix[ia] = (astr.crpix[0]-1.0)*ratio[ia] + 1.0      else $
     	 newastr.crpix[ia] = (astr.crpix[0]-0.5)*ratio[ia] + 0.5

		; TODO redo this to check for PC version too
		if (numparams EQ 2) then begin
			 ;CDn_m Matrix format
			ratio_matrix = rebin(transpose(ratio),3,3)
			newastr.cd = astr.cd * ratio_matrix
		endif else begin    
			; CDELT format
		  	newastr.cdelt = astr.cdelt * ratio
		endelse
	endfor 
	newheader = cubestruct.header
	putast3,  newheader, newastr

	

	return, {cube: newcube, header: newheader, astr: newastr}
	

end