pro dt_factor,dp_sync,roadmap,dt_fac,type,stime,etime
;+
;
;  Routine to calculate the deadtime correction factor by which to
;  multiply the spectral intensities.
;
;  Assumes the correction factor is LIM/ACC (type = 1) or SUM(LIM)/SUM(ACC)
;  (type = 2) where the summation is taken over channels 1 & 2 or 3 & 4. The
;  counts are corrected for their own deadtime before use. 
;
;  Must have run test_rd first and requested reading of dp_sync data
;  or alternatively used the routines:
;
;        IDL> rd_roadmap,infil,roadmap
;        IDL> rd_bda_dp,infil,dp_sync
;
;             where infil = data file name
;
;  Input:
;           dp_sync     usual structure
;           roadmap     usual structure
;           [type]      define form of correction (1 or 2) default is 1 if not
;                       supplied.
;           [stime, etime]  optional hours UT within which factors are found
;                           if not given then whole file is processed (can
;                           take a long time)
;
;  Output:
;           dt_fac(4,n)     deadtime factors by which to multiply spectral
;                           counts. The dimension n is kept the same as
;                           roadmap for ease of application even if not all
;                           elements were requested.
;
;  Use:
;
;        IDL> dt_factor, dp_sync, roadmap, dt_fac [type, stime, etime]
;
;
;  CDP March 92   Written
;      Jan   93   Set factor on error to 1.0 not 0.0
;
;-

;
;  was method supplied?
;
if n_params() lt 4 then type = 1
if type gt 2 or type lt 1 then begin
   print, 'Invalid method requested - only 1 or 2 allowed.'
   return
endif
;
;  set default limited events counter deadtime
;
lim_microsecs = 3.7
;
;  were times supplied? if not use all file, if yes then find indices
;  in roadmap arrays
;
if n_params() lt 6 then begin
   n1 = 0
   n2 = n_elements(roadmap.time) - 1
endif else begin
   ttime = roadmap.time/3600./1000. + (roadmap.day - roadmap(0).day)*24.
   ss = where( (ttime ge stime) and (ttime le etime) )
   if n_elements(ss) ge 2 then begin
      n1 = ss(0)
      n2 = ss(n_elements(ss)-1)
   endif else begin
      print,'No datasets found in that time range, returning.'
      print,'PS. time range in this file is ',ttime(0),' to ',$
                   ttime(n_elements(roadmap.time)-1),$
                   ' Hrs UT',format='(a,f6.3,a,f6.3,a)'
      return
   endelse
endelse
;
;  and the same in the dp_sync arrays
;
if n_params() lt 6 then begin
   d1 = 0
   d2 = n_elements(dp_sync.time) - 1
endif else begin
   ttime = dp_sync.time +(dp_sync.day - dp_sync(0).day)*86400L*1000L
   ttime = ttime/3600./1000. 
;
; compensate for different start days
;
   if roadmap(0).day lt dp_sync(0).day then ttime = ttime + 24.0
   if roadmap(0).day gt dp_sync(0).day then ttime = ttime - 24.0

   ss = where( (ttime ge stime) and (ttime le etime) )
   d1 = ss(0)
   d2 = ss(n_elements(ss)-1)
endelse

;
;  use chekov to correct overflows in LIM and ACC counts in all channels
;  and form mean of both counts
;
print,'Checking for counter overflows...'
chekov,3,2,dp_sync,lim4,dum
lim4 = (lim4 + dum)/2
chekov,3,3,dp_sync,acc4,dum
acc4 = (acc4 + dum)/2
chekov,2,2,dp_sync,lim3,dum
lim3 = (lim3 + dum)/2
chekov,2,3,dp_sync,acc3,dum
acc3 = (acc3 + dum)/2
chekov,1,2,dp_sync,lim2,dum
lim2 = (lim2 + dum)/2
chekov,1,3,dp_sync,acc2,dum
acc2 = (acc2 + dum)/2
chekov,0,2,dp_sync,lim1,dum
lim1 = (lim1 + dum)/2
chekov,0,3,dp_sync,acc1,dum
acc1 = (acc1 + dum)/2

;
;  correct for their own deadtime, using only the data requested to save time
;
print,'Correcting LIMited counts for counter deadtime...'
dt_corr,lim_microsecs,lim4(d1:d2),dum
lim4_dt = fltarr(n_elements(lim4))
lim4_dt(d1:d2) = dum

dt_corr,lim_microsecs,lim3(d1:d2),dum
lim3_dt = fltarr(n_elements(lim3))
lim3_dt(d1:d2) = dum

dt_corr,lim_microsecs,lim2(d1:d2),dum
lim2_dt = fltarr(n_elements(lim2))
lim2_dt(d1:d2) = dum

dt_corr,lim_microsecs,lim1(d1:d2),dum
lim1_dt = fltarr(n_elements(lim1))
lim1_dt(d1:d2) = dum

;
;  form the search index array
;
rm_idx = indgen(n2-n1+1) + n1

;
;  call routine to correlate roadmap and dp_sync times
;
dp_rmap_time,dp_sync,roadmap,rm_idx,dp_idx
;
;  create the output factor array
;
dt_fac = fltarr(4,n_elements(roadmap.time))
;
;  set all to a safe value of 1.0   *** change ***
;
dt_fac(*,*) = 1.0

;
;  average the counts and form the ratios
;
case type of
   1: print,'Forming type 1 deadtime correction ratios in dt_fac(4,*)...'
   2: print,'Forming type 2 deadtime correction ratios in dt_fac(4,*)...'
endcase
for i=n1,n2 do begin
   if dp_idx(i,0) ge 0 and dp_idx(i,1) ge 0 then begin
      case type of
      1: begin
         dt_fac(3,i) = total(lim4_dt(dp_idx(i,0):dp_idx(i,1)))/$
                       total(acc4(dp_idx(i,0):dp_idx(i,1)))
         dt_fac(2,i) = total(lim3_dt(dp_idx(i,0):dp_idx(i,1)))/$
                       total(acc3(dp_idx(i,0):dp_idx(i,1)))
         dt_fac(1,i) = total(lim2_dt(dp_idx(i,0):dp_idx(i,1)))/$
                       total(acc2(dp_idx(i,0):dp_idx(i,1)))
         dt_fac(0,i) = total(lim1_dt(dp_idx(i,0):dp_idx(i,1)))/$
                       total(acc1(dp_idx(i,0):dp_idx(i,1)))
         end
      2: begin
         dt_fac(3,i) = (total(lim4_dt(dp_idx(i,0):dp_idx(i,1)))+$
                       total(lim3_dt(dp_idx(i,0):dp_idx(i,1))))/$
                       (total(acc4(dp_idx(i,0):dp_idx(i,1)))+$
                       total(acc3(dp_idx(i,0):dp_idx(i,1)))) 
         dt_fac(2,i) = dt_fac(3,i)

         dt_fac(1,i) = (total(lim2_dt(dp_idx(i,0):dp_idx(i,1)))+$
                       total(lim1_dt(dp_idx(i,0):dp_idx(i,1))))/$
                       (total(acc2(dp_idx(i,0):dp_idx(i,1)))+$
                       total(acc1(dp_idx(i,0):dp_idx(i,1)))) 
         dt_fac(0,i) = dt_fac(1,i)
         end
      endcase
   endif else begin
      dt_fac(*,i) = 1.
   endelse
endfor
end