;+
;
; NAME: Enorm_drm
;
; PURPOSE: Scale the DRM output rows from energy units to units of
;	   the peak of the DRM row.  This is an intermediate routine
;	   for scaling scattering matrices to the edges of the direct
;	   matrices
;
; CATEGORY:
;	Spectral Analysis
;
; CALLING SEQUENCE:
;       DRM = $
;	 ENORM_DRM( IN_DRM=IN_DRM, IN_E_OUT=IN_E_OUT, ENORM=ENORM,$
;		ESCALE=ESCALE, USESPLINE=USESPLINE )
;
; CALLS TO:
;	INTERPOL, SPLINE
;
; INPUTS:
;       IN_DRM 	 - Response matrix in cnts/energy
;	IN_E_OUT - Output energy edges for IN_DRM,  2 x N, Lo_edge, Hi_edge
;	USESPLINE   - If set, use spline interpolation.
;
; OUTPUTS:
;	ENORM    - Normalized energy scale for output rows
;	ESCALE   - Scaling energy used in each row
; PROCEDURE:
;	Uses INTERPOL, SPLINE, and GAUSSFIT to find the energy at the peak
;	of an output row, then interpolates the row into channels in units of
;	the peak energy.
;
; MODIFICATION HISTORY:
;	Written 19-Aug-1993, RAS
;-

function enorm_drm, in_drm=in_drm, in_e_out=in_e_out, escale=escale, $
	enorm=enorm, usespline=usespline



;Interpolate to geometric average energies

in_em  = (sqrt( in_e_out(0,*)  *  in_e_out(1,*) ))(*)
enorm  = (findgen(30)+.5) /10. ;0-3, normalized output energies
pre_enorm = -4+findgen(5)
enorm1 = [pre_enorm, enorm]
escale = in_em*0.

ndrm = fltarr( n_elements(enorm), n_elements(in_drm(0,*)) )
for i=0,n_elements( ndrm(0,*) ) -1 do begin 
	evec = in_em/in_em(i)

	if keyword_set(usespline) then $
	  test1 = spline([-1,0, evec, max(evec)+1.+findgen(2)],$
			 [0,0,in_drm(*,i),0,0], enorm) $
	else $

	  test1 = interpol([0,in_drm(*,i)], [0, in_em/in_em(i)], enorm) 

	test2 = gaussfit( enorm1, [fltarr(5),test1]>0, apar)
	escale(i) = apar(1)*in_em(i)

	if keyword_set(usespline) then $
	  ndrm(*,i) = spline([-1,0, evec,$
			max(evec)+1.+findgen(2)]/apar(1),$
				[0,0,in_drm(*,i),0,0], enorm) >0.0 $

	else $

	  ndrm(*,i) = interpol([0,in_drm(*,i)], $
			[0, in_em/(in_em(i)*apar(1))], enorm) >0 

endfor

return, ndrm
end
;----------------------------------------------------------------------------

pro add_scatter_drm, ddrm=ddrm, sdrm=sdrm, new_drm=new_drm, shead=shead, $
	dhead = dhead, nsmooth=nsmooth
;+  
;
; NAME: 
;	ADD_SCATTER_DRM
;
; PURPOSE:
;	Take the scatter matrix, binned at different resolution from the
;	direct matrix, rescale it to the input and output channels of the
;	direct matrix, and add the direct and new scatter matrix into a single
;	detector response matrix.
;
; CATEGORY:
;	BATSE, Spectroscopy
;
; CALLING SEQUENCE:
;	add_scatter_drm, ddrm=ddrm, sdrm=sdrm, new_drm=new_drm, shead=shead, $
;	dhead = dhead [nsmooth=nsmooth]
;
; INPUTS:
;       ddrm - The direct response matrix, see Read_bdas_drm.pro, output
;		rows in counts/cm2/keV
;	sdrm - Another response matrix to add with different input and output
;		binning, usually the atmospheric scattering matrix which 
;		takes too long to produce at the full resolution.
;
;	The following two input structures also come out of 
;	READ_BDAS_DRM.PRO, but they only require the tags described:
;	dhead - Structure with channel binning energy edges in
;		dhead.nobins_in - number of input bins
;		dhead.nobins_out - number of output bins
;		dhead.e_in       - 2 x n input bin edges (keV), lo and hi
;		dhead.e_out      - 2 x n output bin edges
;	shead  - Structure like dhead for scattering matrix
;				
; OPTIONAL INPUTS:
;	nsmooth - average resultant scatter matrix over this number of rows.
;		  The average is taken over the total effective area, and used
;		  to scale the rows.  This removes the scalloped effect of
;		  the interpolations.
;
; OUTPUTS:
;	new_drm - Summed detector response matrix with edges of ddrm
;
; PROCEDURE:
;	Input rows of SDRM are interpolated prior to align them
;	along E_out(*)/e_in(i).  Then the rows are interpolated along
;	the column energies to the inputs of DDRM.  Then finally
;	interpolated to the e_out(*)/e_in(i) values of the edges of DDRM.
;	Once interpolated, then the new SDRM is added to DDRM.
;
; MODIFICATION HISTORY:
;	ras, 18 Aug 93
;-

 
e_in_s=shead.e_in(*,0:shead.nobins_in-1)
e_out_s=shead.e_out(*,0:shead.nobins_out-1)

em_in_s = ( e_in_s(0,*) * e_in_s(1,*) )^.5
em_out_s = ( e_out_s(0,*) * e_out_s(1,*))^.5

e_in_d=dhead.e_in(*,0:dhead.nobins_in-1)
e_out_d=dhead.e_out(*,0:dhead.nobins_out-1)

em_in_d = ( e_in_d(0,*) * e_in_d(1,*) )^.5
em_out_d =( e_out_d(0,*) * e_out_d(1,*))^.5


;change rows of scattering matrix to normalized energies
sdrm_0 = enorm_drm( in_drm=sdrm, in_e_out=e_out_s, $
		usespline=usespline, enorm=enorm, escale=escale) ;
min_sdrm = min(sdrm_0(where( sdrm_0 gt 0)))
sdrm_0 = sdrm_0 > (min_sdrm/1e4)

;interpolate along normalized output energies using input energies
;log/log interpolate along normalized output energies using input energies

sdrm_1 = transpose(x_eout_drm( in_drm= transpose(sdrm_0), $
        	in_e_out=e_in_s, out_e_out = e_in_d,$
		usespline=usespline, /uselog)) 

min_sdrm = min(sdrm_1(where( sdrm_1 gt 0)))
sdrm_1 = (sdrm_1 > (min_sdrm/1e4)) > 1e-25

new_escale = exp(interpol( alog(escale), alog(em_in_s), alog(em_in_d)))

sdrm_2 = ddrm * 0.0 ;

for i=0,n_elements(ddrm(0,*))-1 do $
	sdrm_2(*,i) = exp(interpol( alog(sdrm_1(*,i)), $
	 	alog(enorm*new_escale(i)), alog(em_out_d)))

sdrm_2 = sdrm_2 >0

checkvar, nsmooth, 1
if nsmooth gt 1 then begin
	w_out_d = e_out_d(1,*)-e_out_d(0,*)
	ww = rebin(w_out_d(*),dhead.nobins_out,dhead.nobins_in)
	eff =rebin(ww*sdrm_2,1,dhead.nobins_in) > 1e-20
	smooth_factor = smooth(eff,nsmooth) / eff
	sdrm_2 = sdrm_2 * rebin(smooth_factor,dhead.nobins_out, $
						dhead.nobins_in)
endif

new_drm = ddrm + sdrm_2

return
end