;+
;
; NAME: 
;       PLOT_OCCS
;
; PURPOSE:
;       This procedure controls the extraction of the source occultation steps from the
;       raw data files and displays the fitted model including the occultation step.
;
; CATEGORY:
;       BATSE, OCCULTS
;
; CALLING SEQUENCE:
;	plot_occs, occout=occout, occtime=occtime, $
;	start_time=start_time, end_time=end_time, daynight=daynight, $
;	nsum=nsum, inrange=inrange, tdur=tdur, ord=ord, step=step, $
;	cont=cont, file = file, det_id=det_id, ntop=ntop, tmin=tmin, cut=cut
;
;
; CALLS:
;	regression, stdev, read_occ_discsp, anytim, read_cont, linecolors, utplot
;	setutbase, clear_utplot, datatype
; INPUTS:
;       occout: a structure with the tags
;	start_time, end_time, daynight, occtime, det_id
;
;	or input directly through keywords
;       start_time: start time for dataset, seconds from 1-jan-1979
;       end_time:    end time for dataset
;       daynight:    if 1 then this occultation goes from spacecraft night to day
;       occtime:     midpoint of occultation extinction in seconds from 1-jan-1979
;       det_id:       detector ID
;       
;       nsum:    sample average for display [4]
;       inrange: yrange for display
;       tdur:      include this many seconds around both sides of step in fit [60.]
;       ord:       order of polynomial model used for background [1]
;       tmin:     exclude this many seconds on each side of occultation step [8.]
;       cut:       exclude points from the fit which exceed this number of standard
;                    deviations in the model fit.
;       priority: specifies database order for discsp, default: ['occ_discsp','.discsp','discsp_msfc']
; OPTIONAL INPUTS:
;	none
;
; OUTPUTS:
;       STEP: fltarr(3, nchans)
;	For each channel the occultation step count rate(/sec) returned from the regression as well
;	as the sigma for the fit. Third value is the average count rate (/sec) at the
;	time of the fit.
;
; OPTIONAL OUTPUTS:
;	none
;
; KEYWORDS:
;	noplot - don't plot
; COMMON BLOCKS:
;	plot_occs
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	Currently limited to CONT, DISCLA, and DISCSP data
;
; PROCEDURE:
;	Data are read channel by channel.  For each channel the data are fit to a 
;	model including a polynomial in time for the background as well as
;	a step in the countrate at the time of the occultation excluding the
;	interval 2*tdur centered on the occultation itself.  Outliers are excluded
;	from the fit above cut.
;
; MODIFICATION HISTORY:
;	Version 1, ras, 6-June-1996
;	Version 2, ras, 23-Sept-1996
;	Version 3, ras, 27-Sept-1996 added bdb files for DISCLA
;-
pro plot_occs, occout=occout,step=step, occtime=occtime, $
	start_time=start_time, end_time=end_time, daynight=daynight, $
	nsum=nsum, inrange=inrange, tdur=tdur, ord=ord, file=file, filpath=filpath, $
	cont=cont, det_id=det_id, ntop=ntop, tmin=tmin, cut=cut, priority=priority, $
	noplot=noplot

if n_elements(occout) eq 1 then begin
	start_t=occout.start_time
	end_t  =occout.end_time
	dn     =occout.daynight
	ut_occ =occout.ut_occ
	det    =occout.det_id
endif

if n_elements(start_time) ge 1 then start_t=anytim(/utime,start_time)
if n_elements(end_time) ge 1 then end_t=anytim(/utime,end_time)
if n_elements(occtime) ge 1 then ut_occ=anytim(/utime,occtime)
if n_elements(daynight) ge 1 then dn = daynight
if n_elements(det_id) ge 1 then det = det_id
checkvar, nsum, 4
checkvar, tdur, 60.		;include this many seconds around both sides of step	
checkvar, ord, 1
checkvar, tmin, 8		;exclude this many seconds around both sides of step	
checkvar, cut, 3


if not keyword_set(noplot) then linecolors

start_t = start_t > anytim( ut_occ,/dat,/utime)

common plot_occs, data_in, data_start, data_end, data_tipe, det_id_com, ut

checkvar, data_start, 0
checkvar, data_end, 0
checkvar, data_tipe, 0
checkvar, det_id_com, -1
if occ_filetype(file) eq 'occ_discsp' then request_id=det else request_id='*'
case 1 of
	datatype(request_id) eq 'STR':  det_check = request_id eq '*' and total(det_id_com) eq 28 
	else: 	case 1 of
	n_elements(det_id_com): det_check = request_id eq det_id_com
	else: det_check = 0
	endcase
endcase


if total([abs(data_start - start_t),abs(data_end - end_t)]) ge 1 or data_tipe ne fcheck(cont,0) $
	or not det_check then doread=1 else doread=0

if not keyword_set(cont) then begin
	if doread then begin
		ut=0
		data_in=0
		read_occ_discsp, start_t, end_t, ut, data_in, filpath=filpath, $
		det_id=det_out, priority=priority, file=file, request_id=request_id, /nohkg
		det_id_com = fcheck(det_out,-1)
		data_start=start_t
		data_end  =end_t
		data_tipe = 0 ;discsp or discla data
	endif 

	endif else begin
	if doread then begin
		read_cont, cont= data_in, seconds=ut, file=file, $
			start_t=start_t, end_t=end_t
		data_start=start_t
		data_end  =end_t
		data_tipe = 1 ;discsp data
		data_in = data_in(1:fcheck(ntop,10),*,*)
	endif 

	endelse
clear_utplot
setutbase, ut_occ
sdata_in = size(data_in)
if sdata_in(0) eq 2 then begin
	data_in = reform(data_in,sdata_in(1),1,sdata_in(2))
	idet = 0
endif else idet=det <( n_elements(data_in(0,*,0))-1)


wnz = where(data_in(0,idet,*) ne 0, nnz) 

 
if nnz ge 10 then begin
x = ut(wnz) - ut_occ

!p.multi =[0,1,n_elements(data_in(*,0,0))] 
step = fltarr(3, n_elements(data_in(*,0,0)) )

w = where( abs(x) le tdur and abs(x) ge tmin, nw)
if nw ge 10 then begin
	x = x(w)
	wm = where( x gt 0, nwm)
	if nwm ge 5 and (nw-nwm) ge 5 then begin
	nchan = n_elements(data_in(*,idet,0))
	if nchan gt 2 then charsize=2 else charsize=1
	xmodel = fltarr(3,nw) + 1.
	if dn then xmodel(0,0:wm(0)-1)=0. else xmodel(0,wm)=0.
	xmodel(1,*) = x(*)
	xmodel(2,*) = x(*)^2
	for ichan=0,nchan-1 do begin

		regression, xmodel(0:ord,*), data_in(ichan,idet,wnz(w))/2.048, x*0.0+1.,a0, coeff, resid, $
		 yfit, sigmafit, /noprint
		deviants = yfit(*) - (data_in(ichan,idet,wnz(w))/2.048)(*)
		stdev = sigma(deviants)
		wuse = where( abs(deviants/stdev) le cut, nwuse)
		if nwuse ge 6 then begin
		regression, xmodel(0:ord,wuse), data_in(ichan,idet,wnz(w(wuse)))/2.048, $
			x(wuse)*0.0+1.,a0, coeff, resid, yfit, sigmafit, /noprint
		step(0:1,ichan) = [coeff(0),sigmafit(0)]
		yrange = avg( data_in(ichan,idet,wnz(w(wuse)))/2.048)
		step(2,ichan) = yrange
		yrange = yrange + [-50,50]*sqrt( (yrange/100. > 1.))

		if not keyword_set(noplot) then begin

		utplot, ut-ut_occ, data_in(ichan,idet,*)/2.048, psym=10, nsum=nsum,$
	 	yrange=fcheck(inrange,yrange), ymin=-1, $
		charsiz=charsize, ymar=[2,2],titl='DET '+strtrim(det,2)
		oplot, dblarr(2), crange('y'), color= 7-2*dn

		oplot, xmodel(1,wuse), yfit, psym=2,  color=5

		endif
		endif
	endfor
endif
endif
endif
end