;+
; Project     : SOHO - CDS     
;                   
; Name        : RASTER_DUR()
;               
; Purpose     : Calculates duration of a CDS raster.
;               
; Explanation : Simple calculations of the CDS raster duration based on
;               exposure time and telemetry data size were not sufficiently
;               accurate.  This routine simulates the CDHS data processing
;               and thus models the on-board raster process.  Although
;               it will always be impossible to predict the duration of a
;               raster precisely, this routine mimics the CDHS behaviour 
;               to give durations accurate to a few seconds.  See CDS
;               Software Note #xx for further details.
;               
; Use         : IDL> dur = raster_dur(ras_id=ras_id, ras_varn=ras_varn.....)
;                              
; Inputs      : ras_id - Fundamental raster ID
;               ras_varn - Raster variation ID
;               
; Opt. Inputs : None
;               
; Outputs     : None
;               
; Opt. Outputs: pktb - array with packet buffer history.
;               expst  - array with exposure intervals
;               numexp - number of exposure in raster
;               
; Keywords    : ulatency - overrides calculated latency
;               buff     - sets initial pkt buffer contents
;               erate    - sets engineering rate (default 8 secs)
;               debug    - gives extra info
;               tpi      - telemetry packet interval (msec)
;               pr       - packet processing rate (pkts/sec)
;               spare    - (output) spare inter-exposure telemetry time
;               maxtick  - maximum number of packet time ticks
;               uexp     - overide number of exposures
;               et       - overide raster exposure time
;               dump     - calculate data dump time only
;
; Calls       : Database access routines
;
; Common      : dur_common  - to minimise database calls
;               
; Restrictions: None
;               
; Side effects: None
;               
; Category    : Engineering, planning
;               
; Prev. Hist. : None
;
; Written     : C D Pike, 8-Nov-97
;               
; Modified    : Modify to replace calc_ras_dur for time2dump.  CDP, 12-Dec-97
;               Adjusted GIS latency.  CDP, 16-Dec-97
;
; Version     : Version 3, 16-Dec-97
;-            

function raster_dur,pktb,expst,ras=ras,var=var,debug=debug,tpi=tpi,$
                 ulatency=ulatency,buff=buff,erate=erate,spare=spare,$
                 pr=pr,mechset=mechset,et=et,dfac=dfac,time2dump=time2dump,$
                 maxtick=maxtick,uexp=uexp,ucomp=ucomp,numexp=numexp,dump=dump

;
;  common to save database access
;
common dur_common, dur_state, mech_settle, mech_step

;
;  initialise
;
time2dump = 0.0

;
;  Time interval in secs. between telemetry packets in normal telemetry
;
if keyword_set(tpi) then begin
   pkt_int = tpi/1000. 
endif else begin
   pkt_int = 0.205
endelse


;
;  Fundamental and variation IDs must be given        
;
if keyword_set(ras) and keyword_set(var) then begin
   if datatype(ras,1) eq 'Integer' and datatype(var,1) eq 'Integer' then begin
      errmsg = ''
      get_raster,ras,var,raster,errmsg=errmsg
      if errmsg ne '' then begin
         print,errmsg
         return,0
      endif
      exp_time   = raster.exptime
      nwin       = raster.n_windows
      compid     = raster.comp_id
      compopt    = raster.comp_opt
      vds_orient = raster.vds_orient
      detector   = raster.detector
      nx         = raster.nx
      ny         = raster.ny
      xstep      = raster.xstep
      ystep      = raster.ystep


      get_datawin,raster.dw_id, dw
      if dw.dw_id lt 0 then begin
         print,'Invalid DW_ID"
         return,0
      endif
   endif 
endif else begin
   if keyword_set(ras) and datatype(ras,1) eq 'Structure' then begin
      exp_time   = ras.raster_v.exptime
      compid     = ras.raster_v.comp_id
      compopt    = ras.raster_v.comp_opt
      vds_orient = ras.raster_v.vds_orient
      detector   = ras.raster_p.detector
      nx         = ras.raster_p.nx
      ny         = ras.raster_p.ny
      xstep      = ras.raster_p.xstep
      ystep      = ras.raster_p.ystep
      dwid       = ras.raster_v.dw_id
      nwin       = ras.data_win.n_windows
  
;
;  no data
;
      if nwin eq 0 then return,0

;
;  set packet interval according to input definition
;
      if not keyword_set(tpi) then begin
         case ras.raster_v.tel_rate of
            'L': pkt_int = 1.249
            'P': pkt_int = 0.7171
            'M': pkt_int = 0.205
            'H': pkt_int = 0.1035
           else: pkt_int = 0.205
         endcase
      endif

;
;  Define the output structure as a whole.
;
      dw = {DW_ID:ras.raster_v.DW_ID,         $
              DETECTOR: ras.raster_p.detector,  $
              DW_DESC: ras.data_win.dw_desc,    $
              W_WIDTH: (0),    $
              W_HEIGHT:(0),  $
              VDS_BACK: ras.data_win.VDS_BACK,  $
              N_FLAGS: (0),    $
              LL_ID: ras.raster_v.LL_ID,        $
              LL_DESC: ras.line_list.LL_DESC,    $
              N_LINES: ras.line_list.N_LINES,    $
              WINS: ras.data_win2(0:nwin-1)}
     endif else begin
      print,'Must define raster and variation IDs or give internal structure'
      return,0
   endelse
endelse

;
;  check for silly exposure time that has dire consequences
;
exp_time = exp_time > (2.0*pkt_int)

;
;  user overide on compid?
;
if n_elements(ucomp) gt 0 then compid = ucomp

;
;  check compression is supported
;


case compid of
   1: compfac = 1.0
   2: compfac = 0.75
   6: begin
         if compopt eq 0 then begin
            compfac = 2.0/dw.wins(0).win_def(2)
         endif else begin
            compfac = 2.0/dw.wins(0).win_def(3)
         endelse         
      end
   8: compfac = 1.0/float(compopt)
else: begin
         bell
         print,'Unsupported compression scheme (id: '+trim(compid)+')'
         print,'Assuming truncate compression:'
         compfac = 0.75
      end
endcase
;
;  override exp time
;
if keyword_set(et) then exp_time = et

; 
;  Rate at which packets are created (pkts/sec) duration seems insensitive 
;  to this
;
if keyword_set(pr) then proc_rate = float(pr) else proc_rate = 50.0

;
;  update data extraction windows
;
if detector eq 'N' then begin
   dw = update_dex(dw)
endif else begin
   dw = update_dex(dw,gset=!def_gset_id)
endelse

;
;  calculate VDS extraction windows
;
if n_elements(dur_state) eq 0 then begin
   get_utc,utc
   get_vds_state,utc2tai(utc),dur_state
endif
vdsew = get_vds_win(rotate(dw.wins.win_def,1),dur_state.vds_read)


;
;  calculate number of pixels to read from CCD, hence readout time
;
vtot = 0L
nw = n_elements(vdsew)/4
for i=0,nw-1 do begin
   vtot = vtot + long(vdsew(i,2))*long(vdsew(i,3))
endfor
;
;  readout time from data by Ron Yurow, converted to timer ticks
;
if detector eq 'N' then begin
   tflush = 845.0
   readout = 215.86+vtot*14./1000.+(512L*1024L-vtot)*1.2/1000.+tflush
   off_chip_ticks = round(readout/(pkt_int*1000.0))  
   if keyword_set(debug) then print,'Readout: '+trim(readout/1000.)+' ticks: '+$
                                     trim(off_chip_ticks)
endif else begin
   off_chip_ticks = 0
endelse
init_chip = off_chip_ticks


;
;  and the total number of pixels to be telemetered
;
dextot = 0L
nw = n_elements(dw.wins.win_def)/4
for i=0,nw-1 do begin
   dextot = dextot + long(dw.wins(i).win_def(2))*long(dw.wins(i).win_def(3))
endfor
if keyword_set(dfac) then dextot = long(dextot * dfac)
num_pix = dextot        

;
; number of exposures to take plus 1
;
max_exp = (nx*ny)+1
numexp = max_exp-1
            
if n_elements(uexp) gt 0 then max_exp = uexp + 1
expst = fltarr(max_exp)


;
;  calculate number of telemetry packets needed, taking into account the 
;  compression method and that exposure header has room for 252 bytes and 
;  normal data packets 276 bytes.
;
pbyte = dextot*compfac*2.-252.
if long(pbyte) mod 276 eq 0 then remain = 0 else remain = 1
pkt = long(pbyte)/276 + remain + 1 
if keyword_set(debug) then begin
     print,'Number of pixels/packets to extract = ',dextot,'  ',pkt
     print,'Exposure time: '+trim(exp_time)+'  Dump time: '+$
                trim(pkt*pkt_int,'(f5.1)')
     print,'Number of exposures (nx,ny): ',nx,ny
     if n_elements(uexp) gt 0 then begin
        print,'Overidden by uexp = ',max_exp-1
     endif
endif

;
;  Pure dump time of telemetry
;
time2dump = pkt*pkt_int


;
;  return if /DUMP set
;
if keyword_set(dump) then return,0.0

;
;  setup related variables used later
;
prod_pkt = fix(pkt)           
init_prod_pkt = fix(pkt)     


;
;  arrays needed to keep track of whats happening.
;  For each tick, note contents of packet buffer 
;  (pktb) and also the ticks when exposures start (expst).
;
if n_elements(maxtick) eq 0 then maxtick = 500000L
pktb = intarr(maxtick)
type = bytarr(maxtick)
;expst = fltarr(nx*ny+1)

;
;  accumulation variables for each type of telemetry packet
;
hkc = 0
eab = 0
dat = 0
md = 0
spare = 0.0

;
;  initialise test example
;

num_exp = 1             ; current number of exposure
integ_gone = 0.0        ; cumulative exposure (secs)
data_in_b1 = 0          ; is there data in B1?
data_in_b2 = 0          ; is there data in B2?
exposing = 0            ; exposure being made?

;
;  debugging option
;
if keyword_set(buff) then pkt_buff = buff else pkt_buff = 0            


;
;  get mechanism latencies
;
;if mech_step eq 0 then begin
   ms = cp_get_entry ('cb5lat', 0) 
   mech_settle = ms.active
   ms = cp_get_entry ('cb5lat', 1)
   mech_step = ms.active
;endif

;
;  user override
;
if keyword_set(mechset) then mech_settle = mechset 

;
;  calculate baseline delay for data transfer plus basic CDHS latency
;  Two cases of data size < or > 64K
;
if detector eq 'N' then begin
   if (nx gt 1) and (xstep gt 0) then begin
      mech_settle_time = mech_settle/1000.
      mech_step_time = mech_step/1000.
   endif else begin
      mech_settle_time = 0.0
      mech_step_time   = 0.0
   endelse 
endif else begin
   if (ny gt 1) and (ystep gt 0) then begin
      mech_settle_time = mech_settle/1000.
      mech_step_time = mech_step/1000.
   endif else begin
      mech_settle_time = 0.0
      mech_step_time   = 0.0
   endelse 
endelse  

if keyword_set(debug) then begin
   print,'Mechanism times: ', mech_settle_time, mech_step_time
   print,'Dextot = ',dextot
endif

;
;  determine latency depending on mode to be used
;
if dextot lt 64L*1024L then begin
   if keyword_set(debug) then print,'Case: data < 64K'
   data_size = 0
   latency = dextot/1000.*0.05792 + 0.8 + mech_settle_time + $
                                          mech_step_time*(xstep/2)

endif else begin
   if keyword_set(debug) then print,'Case: data > 64K'
   data_size = 1
   if compid eq 6 then begin
      latency = dextot/1000.*0.08 + 0.75 + mech_settle_time + $
                                          mech_step_time*(xstep/2)
   endif else begin
      latency = 1.0 + mech_settle_time + mech_step_time*(xstep/2)
   endelse
endelse

;
;  adjust for GIS
;
if detector eq 'G' then begin
;   latency = 0.6 + mech_settle_time + mech_step_time*ystep
   latency = 1.1 + mech_settle_time + mech_step_time*ystep
   if mech_settle_time eq 0 then latency = 2.0
endif

if keyword_set(debug) then begin
     if keyword_set(var) then begin
        print,'Raster variation: '+trim(var)+' '+'Latency calculated as '+$
               trim(latency,'(f5.2)')
     endif
endif
;
;  debug option to override calculated value
;
if keyword_set(ulatency) then latency = ulatency

;
;  Latency is applied as if it were a delay in transferring the data
;
init_b1_b2 = round(latency/pkt_int+0.5)
rlat = latency/pkt_int
errlat = round((init_b1_b2 - rlat)*100.)
totknock = 0
b1_b2_ticks = init_b1_b2
norm_b1_b2_ticks = b1_b2_ticks

;
;  save latency as GIS will adjust it
;
normal_latency = latency

;
; time to process B2 data to pkts remember that pkts
; is already the compressed size of data so blow it up again
;
b2_comp_ticks = round((pkt/proc_rate/pkt_int/compfac)+0.5)          
init_b2_comp = b2_comp_ticks

;
; engineering data rate
;
if not keyword_set(erate) then eng_rate = 8. else eng_rate = erate           

;
;  hence calculate interval between EngA/B packets (and assume HK comes
;  every 15 secs)
;
eng_ab = round(eng_rate/pkt_int)
hk = round(15./pkt_int)

;
;  introduce some randomness in when Eng/HK will occur in the simulation
;
ran_eng = randomu(seed,1)*eng_ab & ran_eng = round(ran_eng(0))
ran_hk = randomu(seed,1)*hk      & ran_hk =  round(ran_hk(0))

;
;  Start main simulation loop, each tick is 0.205 (medium rate) secs so 
;  max duration is  11 hrs
;
for tick = 0L,maxtick-2 do begin   

;
;  clock up integration time, it is only reset when an exposure start is
;  allowed
;
   integ_gone = integ_gone + pkt_int   
   if (integ_gone lt exp_time) then begin
      if not exposing then expst(num_exp) = tick
      exposing = 1 
      off_chip_ticks = init_chip     
   endif else begin
      exposing = 0
   endelse
 
;
;  if not exposing then start the readout
;  
   if not exposing then off_chip_ticks = (off_chip_ticks - 1) > 0
   if (not exposing) and (off_chip_ticks eq 0) and $
      (num_exp lt max_exp)  then   data_in_b1 = 1

   if data_in_b1 and (not data_in_b2) then begin
      b1_b2_ticks = (b1_b2_ticks - 1) > 0
      if b1_b2_ticks eq 0 then begin
         data_in_b2 = 1
         b1_b2_ticks = init_b1_b2
         num_exp = num_exp + 1
         num_win = -1
;
;  if < 64 K case then OK to trigger exposure start here because B1 is empty
;  as data are transferred to B2. If it's GIS then check whether this is a 
;  flyback exposure in which case add extra latency.
;
         if data_size eq 0 then begin
            if num_exp lt max_exp then begin
               integ_gone = 0.0 
               md = 0
               if detector eq 'G' then begin
                  if num_exp gt 1 and (num_exp mod ny eq 0) then begin
                     latency = normal_latency + $
                               (ny-1)*ystep*mech_step_time + $
                               2.0*mech_settle_time
                     init_b1_b2 = round(latency/pkt_int+0.5)
                     b1_b2_ticks = init_b1_b2
                  endif else begin
                     latency = normal_latency
                     knock =  randomu(seed,1)*100. 
                     knock = round(knock(0))
                     if knock lt errlat then knock = 1 else knock = 0
                     b1_b2_ticks = norm_b1_b2_ticks - knock
                     if knock eq 1 then totknock = totknock + 1
                  endelse
               endif else begin
                  knock =  randomu(seed,1)*100. 
                  knock = round(knock(0))
                  if knock lt errlat then knock = 1 else knock = 0
                  b1_b2_ticks = norm_b1_b2_ticks - knock
                  if knock eq 1 then totknock = totknock + 1
               endelse
            endif else begin
               num_exp = max_exp
            endelse
            data_in_b1 = 0
         endif
      endif
   endif

   if data_in_b2 then begin

;
;  is there room in the packet buffer?
;
     if pkt_buff lt 500 then begin
        prod_pkt = (prod_pkt - fix(proc_rate*pkt_int) )

;
;  case of not enough data to xfer
;
        if prod_pkt lt 0 then begin
           xfer = fix(proc_rate*pkt_int) + prod_pkt
           prod_pkt = 0
        endif else begin
           xfer =  fix(proc_rate*pkt_int)
        endelse
        pkt_buff = pkt_buff + xfer
;
;  what if buffer overloaded?
;
        excess = pkt_buff - 500
        if excess gt 0 then begin
           pkt_buff = pkt_buff - excess
           prod_pkt = prod_pkt + excess
        endif
     endif


;
;  if all packets have been processed (ie have cleared the appropriate
;  buffer) then we are free to start another exposure
;
     if prod_pkt le 0 then begin
;
;  if > 64 K case then trigger of exposure start (B1 empty) had to wait 
;  till here
;
;        if data_size eq 1 or detector eq 'G' then begin
        if data_size eq 1  then begin
           if num_exp lt max_exp then begin
              integ_gone = 0.0
              md = 0 
              knock =  randomu(seed,1)*100. 
              knock = round(knock(0))
              if knock lt errlat then knock = 1 else knock = 0
              b1_b2_ticks = norm_b1_b2_ticks - knock
           endif else begin
              num_exp = max_exp
           endelse
           data_in_b1 = 0
        endif
        data_in_b2 = 0
        if (num_exp lt max_exp) then prod_pkt = init_prod_pkt
     endif
   endif


;
;  always decrement pkt buffer, but if it was time for Engineering
;  or HK then replace prod_pkt because data wasn't sent - also allow
;  for case when EngA/B and HK coincide (double hits)
;
;   if (round(tick*double(pkt_int)*1000L) mod long(pkt_int*1000)) eq 0 then begin
      pkt_buff = (pkt_buff - 1) > 0
;   endif 

   if ((tick+ran_eng) mod eng_ab) le 1 or $
               ((tick+ran_hk) mod hk) eq 0 then begin
;
;  double hits
;
      if ((tick+ran_eng) mod eng_ab) le 1 and $
                ((tick+ran_hk) mod hk) eq 0 then begin
         increment = 2
      endif else begin
         increment = 1
      endelse
      pkt_buff = pkt_buff + increment
   endif

;
;  if the raster is finished (here defined as all data sent) then jump out
;  of the simulation
;
   if (pkt_buff eq 0) and (num_exp ge max_exp) then goto, jump2

;
;  Collect the last tick for which exposing = 1, this then is the time of 
;  the end of the last exposure, which is what FIND_DURATION records.
;  This should be changed to a better definition - ie the time when a new
;  exposure can kick off = time a new raster can start.
;
   if exposing then last_exp = tick

;
;  accumulate statistics
;
   if tick lt maxtick then begin
      pktb(tick) = pkt_buff

;
;  and allocate a type to the current packet and gather stats but only 
;  for an exposure near the end.
;
      if ((tick+ran_eng) mod eng_ab) le 1 or $
                  ((tick+ran_hk) mod hk) eq 0 then begin
         if ((tick+ran_eng) mod eng_ab) le 1 then begin
            type(tick) = 10 
            if num_exp eq max_exp-3 then eab = eab + 1
         endif else begin
            type(tick) = 12
            if num_exp eq max_exp-2 then hkc = hkc + 1
         endelse
      endif else begin
         if pkt_buff gt 0 then begin
            type(tick) = 5
            if num_exp eq max_exp-3 then dat = dat + 1
         endif else begin
            type(tick) = 1
;            if num_exp eq max_exp-3 then md = md + 1
            if num_exp lt max_exp then md = md + 1
         endelse
      endelse
   endif

endfor            ; end of main simulation loop

;
;  jump to here when raster finished
;
jump2:


if keyword_set(debug) then begin
   print,'Packets: Data, eng, hk md ',dat,eab,hkc,md
   print,'Duration of raster (FIND_DURATION)) = '+trim(round(last_exp*pkt_int))+ ' secs'
   print,'Duration of raster (finished telemetry) = '+trim(round(tick*pkt_int))+ ' secs'
   print,' '
   print,' '
   print,' '
endif

;
;  spare telemetry time
;
spare = md*pkt_int

;
;  trim the output arrays
if tick eq maxtick-1 then begin
   bell,3
   print,'Calculation incomplete, maximum numbers of cycles reached'
endif   
pktb = pktb(0:tick)
type = type(0:tick)

if n_elements(expst) gt 2 then begin
   expst = expst - shift(expst,1)
   expst = expst(2:*)*pkt_int
endif else expst = 0.0


;
;  set mech_step to zero if want database access next time
;
if keyword_set(mechset) then begin
   mech_step = 0
endif


;
;  return with the answer
;
return,round(last_exp*pkt_int)


end