;+
; NAME: cube_collapse 
; PURPOSE:
;	Collapse a cube along the wavelength axis.
;
; INPUTS:
; KEYWORDS:
; 	cropmult=		crop the collapsed image to multiples of this many pixels
; 					(for consistency with what cube_rebin does for odd axis
; 					sizes)
; OUTPUTS:
;
; HISTORY:
; 	Began 2007-07-16 17:40:19 by Marshall Perrin 
;-


FUNCTION cube_collapse, cubestruct, zrange=zrange, silent=silent, cropmult=cropmult, $
		median=median

	waveaxis = cube_getwaveaxis(cubestruct)

	if ~(keyword_set(silent)) then print, "Collapsing along wavelength axis "+strc(waveaxis)
	if keyword_set(median) then print, "Collapsing using MEDIAN"
	sz = size(cubestruct.cube)

	if keyword_set(zrange) then begin
		if ~(keyword_set(silent)) then message,/info, "Using Z range "+printcoo(zrange)
		subscripts = fltarr(3,2)
		subscripts[*,1] = sz[1:3]-1
		subscripts[waveaxis,*] = zrange
		to_collapse = cubestruct.cube[subscripts[0,0]:subscripts[0,1],subscripts[1,0]:subscripts[1,1],$
		                subscripts[2,0]:subscripts[2,1]]
		if keyword_set(median) then collapsed = median(to_collapse, dimen=waveaxis+1) else $
			collapsed = total(to_collapse, waveaxis+1) / (zrange[1]-zrange[0]+1)
	endif else begin
		if keyword_set(median) then collapsed = median(cubestruct.cube, dimen=waveaxis+1) else $
			collapsed = total(cubestruct.cube,waveaxis+1) / sz[waveaxis+1]
	endelse


	if keyword_set(cropmult) then begin
		if ~(keyword_set(silent)) then message,/info, "cropping image size to multiples of "+strc(cropmult)
		multsz = mround(sz,cropmult)
		wa = where(indgen(n_elements(cubestruct.astr.naxis)) ne waveaxis)
		;multsz[waveaxis+1] = sz[waveaxis+1]
		collapsed = collapsed[0:multsz[wa[0]+1]-1, 0:multsz[wa[1]+1]-1]

	endif

	; TODO update astrometry info!
	newhead = cubestruct.header
	astr = cubestruct.astr
	wa = where(indgen(n_elements(astr.naxis)) ne waveaxis)
	newnaxes = n_elements(astr.naxis)-1

	; crop out the right part of the CD matrix. I can't figure out any
	; clever way to do this sans for loop. ??
	newcd = dblarr(newnaxes, newnaxes)
	for ia=0L,newnaxes-1 do begin
		newcd[ia,*] = astr.cd[wa[ia], wa]
	endfor 

    NEWASTR = {NAXIS:astr.naxis[wa], CD: astr.cd, CDELT: astr.cdelt[wa], $
                CRPIX: astr.crpix[wa], CRVAL:astr.crval[wa], $
				CUNIT: astr.cunit[wa], $
                CTYPE: string(astr.ctype[wa]), LONGPOLE: double( astr.longpole),  $
                LATPOLE: double(astr.latpole), PV2: astr.pv2 }
	newheader = cubestruct.header
	putast3,newheader,newastr


	return, {image: collapsed, header: newheader, astr: newastr}
	

end