;+
; NAME:  cubedisp_gridspec
; PURPOSE:
; 		Overplot a grid of spectra on top of a
;
; INPUTS:
; KEYWORDS:
; OUTPUTS:
;
; HISTORY:
; 	Began 2007-07-16 17:35:02 by Marshall Perrin 
;-


PRO cubedisp_gridspec, cubestruct, rebin=rebin, zrange=zrange, _extra=_extra, $
		waverange = waverange,plotrange=plotrange,stop=stop, nogrid=nogrid

	sz = size(cubestruct.cube)
	waveaxis = cube_getwaveaxis(cubestruct)
	if ~(keyword_set(rebin)) then rebin=5



	if keyword_set(waverange) then begin
		xyz2adl, 0, 0, indgen(sz[waveaxis+1]), cubestruct.astr, ra, dec, wavelen
		i0 = find_closest(wavelen, waverange[0])
		i1 = find_closest(wavelen, waverange[1])
		zrange = [i0,i1]
		if ~(keyword_set(silent)) then message,/info, "Using pixels "+printcoo(zrange)+" for wavelenths "+printcoo(waverange)

	endif


	;--- display the 2D image ----
	if ~(keyword_set(silent)) then message,/info, "Collapsing"
	cont = cube_collapse(cubestruct, zrange=zrange, cropmult=rebin)

	imdisp, alogscale(cont.image, 0, 1e-16), /axis, out_pos=op ,/erase ,xticklen=1e-20,yticklen=1e-20
	if ~(keyword_set(silent)) then message,/info, "Rebinning"

	if keyword_set(nogrid) then return


	; --- rebin the 3D cube ---
	cubesmall = cube_rebin( cubestruct, rebin=rebin, zrange=zrange)

	szs = size(cubesmall.cube)

	wa = where(indgen(n_elements(cubestruct.astr.naxis)) ne waveaxis)
	nx = szs[wa[0]+1]
	ny = szs[wa[1]+1]

	if (keyword_set(stop)) then stop

	;--- display the 3D cube ---
	basex = op[0]
	basey = op[1]
	widthx = op[2]-op[0]
	widthy = op[3]-op[1]

	tmpcube = transpose(cubesmall.cube, [waveaxis, wa[0], wa[1]]) ; put wavelen axis first
	if ~(keyword_set(plotrange)) then plotrange = [0,1e-16]
	for ix=0L,nx-1 do begin
	for iy=0L,ny-1 do begin
		position = [float(ix)/nx*widthx + basex, float(iy)/ny*widthy + basey, $
					float(ix+1)/nx*widthx + basex, float(iy+1)/ny*widthy + basey]
		plot, /noerase, tmpcube[*,ix,iy],position=position, xstyle=5, ystyle=5, $
			yrange=plotrange,_extra=_extra
		if ix ne nx-1 and iy eq 0 then plots, [float(ix+1)/nx*widthx + basex,float(ix+1)/nx*widthx + basex], op[[1,3]],_extra=_extra,/normal
		if iy ne ny-1 and ix eq 0 then plots, op[[0,2]], [float(iy+1)/ny*widthy + basey,float(iy+1)/ny*widthy + basey],_extra=_extra,/normal
	endfor 
	endfor 

	if (keyword_set(stop)) then stop

end