;+
; NAME: sotsp_rasterize
;
; PURPOSE: Takes a list of slit positions, and their observation times, from
;          SOT-SP level 2 data files and returns a mapping file that allows
;          the user to reconstruct the scan (single-scan case) or raster
;          (repeating scan case). Output is a vector or array containing
;          indices that map to the input list of slit positions, accounting
;          for duplicates and missing columns. The routine first looks at the
;          observation times associated with each slit position and screens
;          out duplicates; those with an unrealistic interval of time (see
;          dupe_thresh variable below) between successive scans are
;          removed. Slit positions that are missing are indicated with a value
;          of -9999 in the output array.
;
; CALLING SEQUENCE:
; slitmap=sotsp_rasterize(scan_info,duplicate_indices=duplicate_indices,error=error,/debug)
;
; INPUTS: scan_info - scan_info structure as returned from read_sotsp
;
; KEYWORD PARAMETERS:
;          duplicate_indices - on output, returns an array of indices of the
;                              input slit position list whose slit positions
;                              are too close in time to the one immediately
;                              preceeding it, and are thus considered
;                              duplicates (if none found, then -1 is returned)
;          error - on output, this keyword returns an integer corresponding to
;                  an error code: 0 = no error 1 = problem with finding slit
;                  position in slit index
;          debug - a switch that, if set, will stop when an error is detected
;
; OUTPUTS: slitmap - a vector or array that indexes slit positions in input
;
; EXAMPLE: If on input,
;
;          scan_info.slitpos = [300, 302, 304, 306, 308, 300, 302, 306, 308],
;
;          which represents a small raster of two scans with 5 slit positions,
;          then this routine will return the array
;
;          slitmap = [ [ 0, 1,     2, 3, 4],
;                      [ 5, 6, -9999, 7, 8] ],
;
;          which is a 2D array that indexes the raster from the 1D input
;          array. Note that slit position 304 is missing in the second scan,
;          and so the value of -9999 is used to indicate that data for this
;          slit position is missing in the corresponding level 2 data arrays.
;
; MODIFICATION HISTORY:
;   2020 Aug 27 - M.DeRosa - created
;   2022 May 24 - M.DeRosa - added error and debug keywords
;
;-

function sotsp_rasterize,scan_info,duplicate_indices=duplicate_indices,error=error,debug=debug

  ;  set parameters
  dupe_thresh = 0.6  ;  in seconds, used for detecting duplicate columns
  error=0
  
  ;  first screen out duplicates slit positions
  slitpos_in=reform(round(scan_info.slitpos))
  slittime=anytim(scan_info.times)
  timediff=slittime(1:-1)-slittime(0:-2)  ;  difference is in seconds
  whvalid=[where(timediff gt dupe_thresh,nvalid),n_elements(slitpos_in)-1]
  nvalid++  ;  accounts for final slit position (last element above)
  duplicate_indices=where(timediff le dupe_thresh)  ;  for output
  slitpos=slitpos_in(whvalid)

  ;  determine unique list of observed slit positions
  uniqslit=slitpos(uniq(slitpos,sort(slitpos)))
  nuniqslit=n_elements(uniqslit)
  minslit=min(uniqslit,max=maxslit)

  ;  infer interval between slit positions, as commanded
  slitdiff=slitpos(1:-1)-slitpos(0:-2)
  slithist=histogram(slitdiff(where(slitdiff ge 0)),omin=omin,omax=omax)
  buffer=max(slithist,whmax)
  slitstep=omin+whmax
  if slitstep eq 0 then slitstep++  ;  set to 1 for sit and stare case
  
  ;  determine list of commanded slit positions, including gaps
  nslitpos=1+round(float(maxslit-minslit)/slitstep)
  slitix=minslit+lindgen(nslitpos)*slitstep

  ;  determine number of repeats
  slitreturns=where(slitdiff le 0,nslitreturns)
  nrepeat=nslitreturns+1
  
  ;  create slitmap
  slitmap=make_array(dim=[nslitpos,nrepeat],/long,value=-9999)
  returncount=0  & prevslitpos=minslit-1  ;  returncount increments each row
  for jj=0,nvalid-1 do begin
    currslitpos=slitpos[jj]
    if currslitpos le prevslitpos then returncount++
    whslit=where(currslitpos eq slitix,nwh)
    if nwh eq 0 then begin  ;  error
      error=1
      if keyword_set(debug) then stop,'  problem with finding slit position in slit index'
    endif else slitmap[whslit[0],returncount]=whvalid[jj]
    prevslitpos=currslitpos
  endfor

  ;  return
  return,slitmap

end