;+
; NAME:  cube_subtractcontinuum
; PURPOSE:
; 		Crude continuum subtraction for IFS datacubes
;
; 		Almost any algorithm will be better than this!
;
; INPUTS:
; KEYWORDS:
;
; 	/mean		use a crude moving average to compute the continuum
; 	/median		use a moving median to compute the continuum
; OUTPUTS:
;
; HISTORY:
; 	Began 2007-07-13 17:19:22 by Marshall Perrin 
;-

PRO cube_subtractcontinuum, cubestruct, subtracted, contcube=contcube, $
	median=median, mean=mean

	if ~(keyword_set( median)) and ~(keyword_set(mean)) then mean=1

	; setup: rearrange cube to have wavelength axis in the first position
	; (fastest)
	sz = size(cubestruct.cube)
	axis_lambda = cube_getwaveaxis(cubestruct)
	
	axis_other = (indgen(3))[ where(indgen(3) ne axis_lambda) ]
	tmpcube = transpose(cubestruct.cube, [axis_lambda, axis_other] )
	sz = size(tmpcube)
	z = sz[1]

	; now create a continuum cube
	if keyword_set(mean) then begin 
		stepsize = 50
		zr = mround(z,stepsize)
		tmpcube = tmpcube[0:zr-1,*,*]
	
		print, "binning down"
		cont = rebin(tmpcube, zr/stepsize,sz[2],sz[3] )
		print, "binning up"
		cont = rebin(cont, zr, sz[2],sz[3])
		print, "subtracting"
			
	endif else if keyword_set(median) then begin 
		stepsize = 50
		cont = fltarr(sz[1],sz[2],sz[3])
		nchunks = ceil(z*1.0/stepsize)
		print, "moving median, boxsize = ",stepsize
		for i=0L,nchunks-1 do begin
			statusline, "chunk "+strc(i)
			tmpim = median(tmpcube[ i*stepsize: ((i+1)*stepsize)<z, *,*], dim=1)
			cont[ i*stepsize: ((i+1)*stepsize)<z, *,*]  =  rebin(tmpim,  (((i+1)*stepsize)<z) -(i*stepsize), sz[2], sz[3])
		endfor 
	end
	
	; Subtract, and re-arrange back into the original order
	tmpcube -= cont
	tmpcube = transpose(tmpcube, [axis_lambda, axis_other])
	
	
	

	subtracted = {cube: tmpcube, header: cubestruct.header, astr: cubestruct.astr}
	if arg_present(contcube) then contcube = {cube: contcube, header: cubestruct.header, astr: cubestruct.astr}

	

end