;+
;
; NAME: 
;	DISCSP_BANDS
;
; PURPOSE:
;	This procedure integrates the DISCSP channels from multiple
;	detectors into a single data stream.
;
; CATEGORY:
;	BATSE, SPEX
;
; CALLING SEQUENCE:
;	DISCSP_BANDS, INPUT, UT, FLUX, EFLUX
;
; CALLS:
;	READ_4_SPEX, CHECKVAR, ARR2STR, DET_COSINES, DISCSP_EDGES, BATSE_SPEC_DRM
;	LOC_FILE, EDGE_PRODUCTS
;
; INPUTS:
;       Input - Flare number from SDAC BATSE archive
;		or a structure with tags:
;		FILENAME - DISCSP data file
;		START_TIME - Accumulation start time (readable by ANYTIM)
;		END_TIME   - Accumulation end time
;		COS_SPEC - 8 direction cosines, use if flare value is 0
;
;
; OUTPUTS:
;	UT   -2xnbins, acc. time interval in seconds since 1 jan 1979
;	FLUX - counts in differential band
;	EFLUX- Poisson uncertainty on counts
;	
; KEYWORDS:
;
; OPTIONAL INPUTS:
;	
;	These values can be entered by the tags in the structure Input or Directly.
;	The values in the structure take precedence if found.
;		FILENAME - DISCSP data file
;		START_TIME - Accumulation start time (readable by ANYTIM)
;		END_TIME   - Accumulation end time
;		COS_SPEC - 8 direction cosines, use if flare value is 0
;	MIN_BAND - optional, minimum energy (keV) difference used in creating differential bands
;	OVERLAP  - Differential bands may overlap edges, useful with MIN_BAND constraint.
;	GETDRM   - Just get the DRM info.
; OPTIONAL OUTPUTS:
;	LTIME - Interval livetimes in seconds, nchan x nbins
;	DELTA_LIGHT - Band width used to normalize DRM channels, fltarr(nchan)
;	AREA - detector geometric area, 127. cm2
;	EDGES- differential PHA energy edges
;	DRM -  detector response matrix, cnts/cm2/keV per photons/cm2/bin
;	ERROR- set if there is an unresolvable problem
;	WCHAN- valid channel indices, nominally all of them
;	DET_IDS-   array of detector ids meeting requirements
;	SCALES- normalization constants used to create differential bands
;	EORD-  band order used
;	TITLE- 'BATSE DISCSP MULTI-BAND '+ARR2STR(STRTRIM(IDS,2),'')
; COMMON BLOCKS:
;	none
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	none
;
; PROCEDURE:
;	This procedure first identifies the datafiles, then reads them using read_4_spex, and
;	then creates a merged drm using build_spec_drm.  The PHA channels are built using 
;	adjacent DISCSP edges and taking the difference between integral channels where
;	the high edge is always greater than 1 MeV.  Normally the adjacent edges lie in
;	different detectors, but the rates and DRM's are scaled to ensure that the DRM column
;	doesn't fall below zero for this new differential band.
;	Flux and DRM for this channel are similarly built and finally the
;	difference DRM channel is normalized into the pseudo-PHA width of the newly created
;	differential band.
;
; MODIFICATION HISTORY:
;	Version 1, RAS, 31-oct-1996
;-
pro discsp_bands, input, ut, flux, eflux, ltime=ltime, delta_light=delta_light, $
	min_band=min_band, overlap=overlap, title=title, bdrm=bdrm, scales=scales,$
	area=area, edges=edges, cos_spec=cos_spec, drm=drm, error=error, wchan=wchan, $
        filename=filename, det_ids=ids, start_time=start_time, end_time=end_time, $
	eord=eord, oedges=edges2


common discsp_bands, drmsave
checkvar, min_band, 0.1 


error = 0

if datatype(input) eq 'STC' then begin
        tags=tag_names(input)
        list=['FILENAME','COS_SPEC', 'START_TIME','END_TIME']
        for i=0,n_elements(list)-1 do if (where(list(i) eq tags))(0) ge 0 then $
        eresult=execute(list(i)+"=flare.(w(0))")
	edge_time = avg( anytim(/sec, [start_time,end_time]))
        flare=0
        endif else begin
        flare = input	
	read_flare, flare, fld
	edge_time=fld.peak_secs
        sflare = strtrim(flare, 2)
	filename = loc_file('*.discsp_'+sflare,path='BATSE_DISCSP')
        endelse
if flare ne 0 then det_cosines,fla=flare,cos_spec=cos_spec, spec_sort=spec_sort,/sort

edges=discsp_edges(edge_time)
;
; Identify the detectors where the lowest edge is below 40 keV and theta is < 70 degrees
;
theta = !radeg*acos(cos_spec)
wuse  = where( edges(0,*) le 40.0 and theta le 70.0, nuse)
if nuse eq 0 then begin
	error=1
	printx, 'Energy ranges and detector cosines out of valid ranges'
	printx, 'Valid:  where( edges(0,*) le 40.0 and theta le 70.0 )
endif
ids= (indgen(8))(wuse)
theta = theta(wuse)
edges = edges(0:1,wuse)
merge_batse_photons, e_in
nph   = n_elements(e_in(0,*))
bdrm = fltarr(2,nuse, nph)
area = 127.
flux = 0
for i=0,nuse-1 do begin
	edge_products,[edges(*,i),1000.],edges_2=eout, width=wout
	batse_spec_drm,edges_out=eout,theta=theta(i),det_id=ids(i),phtn_e=e_in,drm=drmi
	bdrm(*,i,*)=drmi * rebin(wout,2,nph) * area
	bdrm(0,i,*)=bdrm(0,i,*)+bdrm(1,i,*)
  endfor 

for i=0,nuse-1 do begin
	read_4_spex, flare=flare,instrument='BATSE', dformat='DISCSP', files=filename,id_style='id',$
	det_id=ids(i), flux=fluxi, eflux=efluxi, ut=uti, units=units, $                     
	ltime=ltimei, edges=edgesi,  wchan=wchani, start_time=start_time, end_time=end_time, $          
	delta_light=delta_lighti                                                                        
	checkvar, nbins, n_elements(uti(0,*))                                                           
	nchans = 2
	                                                                               
;                                                                                               
;Convert from cnts/cm2/s/kev to cnts/s/kev                                                      
;                                                                                               
	if not keyword_set(flux2) then begin                                                             
	    flux2=fltarr( 2,nuse,nbins)                                                                  
    	    eflux2 = flux2                                                                                
	    ut  =uti                                                                                    
	    ltime2 = flux2                                                                                
    	    endif                                                                                       
                                                                                                
	 flux2(*,i,*)     = fluxi * rebin(delta_lighti, nchans,nbins) * ltimei *area 
	eflux2(*,i,*)    = efluxi * rebin(delta_lighti, nchans,nbins) * ltimei *area
	ltime2(*,i,*)    = ltimei
endfor

flux2(0,*,*) = flux2(0,*,*) + flux2(1,*,*)
;
; Reformat the channel arrays and add a zero valued channel to the flux and BDRM
; which represents a channel that ranges from 1000. keV to 1000. keV
;
edges2 = [reform(edges, nuse*2), 1000.]
edges =  fltarr(2,nuse*2)
drm   = fltarr(nuse*2,nph)
bdrm  = [reform( bdrm, nuse*2, nph), fltarr(1,nph)]
flux2  = [reform( flux2, nuse*2,nbins), fltarr(1,nbins)]
flux  = flux2(0:nuse*2-1,*)
ltime2 = [reform( ltime2, nuse*2,nbins), fltarr(1,nbins)]
ltime  = ltime2(0:nuse*2-1,*)
eflux = flux
scales=fltarr(nuse*2)
;
; Create differential channels using consecutive energies
; Adjust the scaling so the difference DRM channel is greater than zero over
; the range from channel low edge to 100 keV

eord = sort(edges2(*))
help, min_band
for i=0,2*nuse-1 do begin
	i1 = (where( edges2(eord) - min_band ge edges2(eord(i))))(0)
        edges(0,i)= [edges2(eord(i)),edges2(eord(i1))]
        wein   = where( e_in(0,*) ge edges(eord(i)) and e_in(0,*) le 100.0, nein)
	scale  = 1
	;Make sure DRM doesn't go negative in the important range
	if nein ge 1 then scale = 1.1*max((f_div(bdrm(eord(i1),*), bdrm(eord(i),*)))(wein)) > 1. 
	scales(i)=scale
	if eord(i) mod 2 eq 0 and eord(i1)-eord(i) eq 1 then begin
	;This is the case where the next channel is just DISCSP2 above DISCSP1
	;for the same detector.
	drm(i,*)   = bdrm(eord(i),*) - bdrm(eord(i1),*)
	flux(i,*) = flux2(eord(i),*) - flux2(eord(i1),*)
	eflux(i,*)= sqrt(flux(i,*))
	endif else begin
	drm(i,*)   = bdrm(eord(i),*) - bdrm(eord(i1),*)/scale
	flux(i,*) = (f_div(flux2(eord(i),*),ltime2(eord(i),*)) - $
	            f_div(flux2(eord(i1),*),ltime2(eord(i1),*))/scale )*ltime2(eord(i),*)
	eflux(i,*)= (sqrt( f_div(flux2(eord(i),*),ltime2(eord(i),*)^2) + $
	                  f_div(flux2(eord(i1),*),ltime2(eord(i1),*)^2)/scale^2))*ltime2(eord(i),*)
	endelse
endfor

delta_light = (edges(1,*)-edges(0,*))(*)
nchan       = n_elements(delta_light)
drm = drm / area / rebin(delta_light, nuse*2, nph)

if not keyword_set(overlap) then begin
	wchan = [0]
	istart=0
	while istart lt (nchan-1) do begin
	wnext = where( edges(0,*) ge edges(1,istart) )

	istart = wnext(0)
	wchan =[wchan, istart]
	endwhile
	if n_elements(wchan) lt n_elements(flux(*,0)) then begin
		drm = drm(wchan,*)
		flux= flux(wchan,*)
		eflux=eflux(wchan,*)
		ltime=ltime(wchan,*)
		edges=edges(*,wchan)
		delta_light=delta_light(wchan)
	endif
endif

wchan = indgen(n_elements(flux(*,0)))
drmsave = {drm:drm, edges_in:e_in, edges_out:edges, theta:theta, area:area, det_ids:ids,$
	delta_light:delta_light}

title = 'BATSE DISCSP MULTI-BAND '+arr2str(strtrim(ids,2),'')

end