;+
;
; NAME: X_EOUT_DRM
;
;
; PURPOSE:
;	Given a Detector Response Matrix, (cnts/keV), calculated on output
;	energy bins, IN_E_OUT, interpolate the output rows to the actual 
;	detector output bins.
;
; CATEGORY:
;	Spectral Analysis
;
; CALLING SEQUENCE:
;       DRM = $
;	 X_EOUT_DRM( IN_DRM=IN_DRM, IN_E_OUT=IN_E_OUT, OUT_E_OUT=OUT_E_OUT, $
;		NSUB_BINS=NSUB_BINS, USESPLINE=USESPLINE )
;
; CALLS TO:
;	INTERPOL
;
; INPUTS:
;       IN_DRM 	 - Response matrix in cnts/energy
;	IN_E_OUT - Output energy edges for IN_DRM,  2 x N, Lo_edge, Hi_edge
;	OUT_E_OUT- Output energy edges for new DRM
;	[NSUB_BINS] - Divide each output bins into N sub-bins, interp, then sum
;                     Default is 1 sub-bin
;	USESPLINE   - If set, use spline interpolation.
;	USELOG      - If set, use log-log interpolation for linear or spline
; PROCEDURE:
;	Uses INTERPOL
;
; MODIFICATION HISTORY:
;	Written 21-May-1993, RAS
;-

function x_eout_drm, in_drm=in_drm, in_e_out=in_e_out, $
	out_e_out=out_e_out, nsub_bins=nsub_bins, usespline=usespline, $
	uselog=uselog

checkvar, nsub_bins, 1
nsub_bins = nsub_bins > 1
n_in = n_elements( in_e_out(0,*) ) ;number of initial output bins

e_o = out_e_out
num = n_elements(e_o(0,*))
new_num = nsub_bins * num
e_o = interpol( [(e_o(0,*))(*), e_o(1,num-1)], new_num + 1)

e_o = transpose( [[e_o(0:new_num-1)],[e_o(1:*)]] )

;Interpolate to geometric average energies

out_em = sqrt( e_o(0,*) * e_o(1,*) )
in_em  = sqrt( in_e_out(0,*)  *  in_e_out(1,*) )
nin = n_elements(in_em) -1
nout= n_elements(out_em) -1


ndrm = fltarr( new_num, n_elements(in_drm(0,*)) )
for i=0,n_elements( ndrm(0,*) ) -1 do begin 
	wi = where( in_drm(*,i) gt 0.0, ni) 
	if ni ge 1 then begin
	  ni = wi(ni-1)+3 < nin
	  wi = indgen(ni)
		wo = where( e_o(0,*) le in_e_out(1,ni-1 ), numo )
		if numo ge 1 then begin
                 	wo = indgen( wo(numo-1) + 3 < nout)
			if keyword_set(usespline) then $
				if not keyword_set(uselog) then $
					ndrm(wo,i) = $
					spline(in_em(wi), in_drm(wi,i), out_em(wo)) $
				else $
					ndrm(wo,i) = $
					exp( spline(alog(in_em(wi) > 1e-15), alog(in_drm(wi,i)),alog(out_em(wo)))) $
			else $
				if not keyword_set(uselog) then $
               				ndrm(wo,i) = interpol( in_drm(wi,i),in_em(wi), out_em(wo)) $
				else $
					ndrm(wo,i) = $
					exp(interpol( alog(in_drm(wi,i) > 1e-15),alog(in_em(wi)),alog(out_em(wo))))
		endif
	endif
endfor

;Sum over the fine output bins
ndrm = rebin( ndrm, num, n_elements(ndrm(0,*))) 

return, ndrm
end