;+
; NAME:
;     LIMCHEK
; PURPOSE:
;     Main routine to analyze LIMPOINT data.  Responds to scans
;     of type PNTCAL (scan code 8) or TEST DATA (scan code 1).
; CATEGORY:
;     OVRO APC DATA ANALYSIS
; CALLING SEQUENCE:
;     limchek[,filename][,hrec | ,after=after][,debug][,/cmdfile]
; INPUTS:
;     filename   the name of the file containing PNTCAL data.  If
;                  omitted, the file DAILY.ARC in directory
;                  !DEFAULTS.WORKDIR is assumed.
;     hrec       the record number of the SOLAR scan header.  Data
;                  will be processed up to the next EOS segment.
;                  If this argument is given, then AFTER keyword
;                  is ignored.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     after      an optional time string of the standard form
;                  [yyyy.ddd ]hh:mm[:ss] after which to start
;                  looking for a valid PNTCAL or TEST DATA scan.
;     cmdfile    a switch that indicates that the routine was called
;                  from a command file, so issue no modal messages, and
;                  print any error messages to <WORKDIR>\ACAL.MSG.
;     debug      a switch that activates some debugging statements
;                  for getting an indication of intermediate steps.
;                  This switch is not compatible with /CMDFILE switch.
; ROUTINES CALLED:
;     openarc, getdata, lasthrec, tl_decode, decode, limfit
; OUTPUTS:
; COMMENTS:
;     There will have to be modifications to this routine when
;     additional antennas are available.
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 20-Dec-1998 by Dale Gary
;     22-Dec-1998  DG
;       Added proper gain/nonlinearity calibration.  Eliminated
;       saturated data.
;-
pro limchek,filename,hrec,after=after,debug=debug,cmdfile=cmdfile
   common freq,frq
   ; If the /cmdfile switch is set, open a "new" LIMPNT.MSG before doing
   ; anything.  If it exists it will be cleared.
   msglun = 0   ; Indicates LIMPNT.MSG does not exist
   IF (keyword_set(cmdfile)) THEN $
                  openw,msglun,/get_lun,!defaults.workdir+'LIMPNT.MSG'

   ; If no filename, assume standard DAILY.ARC file
   IF (n_elements(filename) EQ 0) THEN $
                  filename=!defaults.workdir+'DAILY.ARC'

   ; Open the file
   lun = openarc(filename,a,nrec)

   IF (n_elements(hrec) EQ 0) THEN BEGIN
      IF (keyword_set(after)) THEN hrec = lasthrec(a,nrec,tstr=after) $
                              ELSE hrec = lasthrec(a,nrec)
   ENDIF

   ; Verify that the specified record has either a PNTCAL or TEST DATA
   ; scan code
   rec = hrec - 1
   data = getdata(rec,a)
   tls = tl_decode(data)
   IF (tls.scancode NE !SCAN.TEST AND tls.scancode NE !SCAN.LIMPOINT) THEN BEGIN
      errstr = 'LIMCHEK: Header found is of wrong type.'
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
      return
   ENDIF

   ; Call routine to set up structures at the start of a new scan.
   ; Returns the record number corresponding to the first data record.
   rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm)
   if (rec eq -1) then begin
      errstr = 'LIMCHEK: Error in file structure or missing segment.'
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
      return
   endif

   nant = header.nant

   ; The PNTCAL scan has NBLK blocks of data, each of duration NOBS seconds.
   ; These parameters are set by the TRAJECTORY.  Currently, the defined
   ; trajectory for LIMPOINT is 83 TRAJECTORY, which has 20 blocks, 15 s per block.
;   hodo = where(*traj.pcmd eq 'HODO   ',nblk)   ; indexes where CMD = 'HODO'
;   if (nblk eq 0) then begin
;      errstr = 'LIMCHEK: Invalid TRAJECTORY for LimPoint'
;      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
;                                ELSE ans = widget_message(errstr,/ERROR)
;      return
;   endif
;   ; NBLK should now be 20
;   tot_time = total((*traj.pdur)(hodo))/2  ; trajectory durations are in 0.5 s units
;   nobs = fix(tot_time/nblk)
   ; Hard-wire the trajectory parameters NBLK = 20, NOBS = 15.  The above code
   ; does this, but is not general enough to be useful.
   nblk= 20
   nobs = 15
   npersec = 1000/cfg.sampintms   ; Number of samples per second
   koff = 2                       ; Number of samples before DOSEQ starts

   ; We allow NAQC seconds to acquire the source, and allow NTOL seconds
   ; on either end for "sloppy timing."
   nacq = 5
   ntol = 2

   ; Length of observing sequence
   lseq = obseq.nx

   nndw = header.ndbw + header.nfgw + 1     ; Number of non-data words
   nws = header.nws                         ; Number of words per sample
   nsr = header.nsr                         ; Number of samples per record
   nsamp = nblk*nobs*npersec                ; Number of samples in entire traj
   bigdata = fltarr(nant,nsamp)             ; Big enough to hold all of the data
   gdata = lonarr(nsamp)                    ; Parallel array of gain data

   ; The number of records to expect should be nrecs=nsamp/header.nsr
   ; (24*48*10/33 = 349).  Note that there should be some samples in the 350th
   ; record, but CPC doesn't send it.
   nrecs = nsamp/header.nsr
;   rec = rec + 1
   datoff = header.rdoff+nndw
   ndat = nant

   ; A late start may mean 1 fewer records, so we must be able to withstand
   ; a short scan. Use ON_IOERROR
   ON_IOERROR,EOD

   ; Read data from the file into a single large array, BIGDATA, of
   ; size (3,24*48*10), which is just a long list of 3 total power channels
   ; Loop over the 349 records
   for i = 0,nrecs-1 do begin
      data = getdata(rec,a)                      ; Read a record
      ; Loop over the 33 samples in each record
      for j = 0, header.nsr-1 do begin
         soff = j*header.nws   ; Additional offset for each sample
         bigdata(*,i*header.nsr+j) = data(datoff+soff:datoff+soff+ndat-1)
         gdata(i*header.nsr+j) = w2l(data,datoff+soff-nndw+1)
      endfor
      rec = rec + 1                      ; Increment record number
   endfor

EOD:
   ON_IOERROR,NULL

   gtable = fltarr(8,nant)

   ; Fill in gain table with nominal attenuator values (10^findgen(8)/2.),
   ; corrected for the attenuation factors from the GCAL
   for i = 0, nant-1 do begin
      attn = (*gparm.pgnp)[i].attn   ; 3-element array of factors
      afac = [  1.0, attn(0), attn(1), attn(0)*attn(1), attn(2), $
                       attn(0)*attn(2), attn(1)*attn(2), attn(0)*attn(1)*attn(2)]
      gtable(*,i) = 10^(findgen(8)/2.)/afac     ; Apply correction from GCAL
   endfor

   ; Set saturated data to NaN
   bad = where(bigdata eq '7fff'x,nbad)
   if (nbad ne 0) then bigdata(bad) = !values.f_nan

   ; Apply gain corrections
   for i = 0, nant-1 do begin
      ; Subtract offsets
      bigdata(i,*) = bigdata(i,*)-(*gparm.poffsets)[i]
      ; Apply nonlinearity correction
      bigdata(i,*) = bigdata(i,*)*(1+(*gparm.pgnp)[i].nlparm*bigdata(i,*)/2047.)
      ; Apply gain (DB setting) corrections
      j = ishft(gdata,-4*i) and '7'x   ; Array of indexes into GTABLE
      bigdata(i,*) = bigdata(i,*)*gtable(j,i)
   endfor
   ; Determine the start and stop indexes into the bigdata array where
   ; good data are expected.
   n = indgen(nblk)
   k1 =  (n*nobs+nacq+ntol)*npersec - 1  +koff  ; Array of end indexes of bad data
   k2 = [0,((n+1)*nobs-ntol)*npersec + 1]+koff  ; Array of start indexes of bad data

   ; Set bad data parts to zero
   for i = 0, nblk-1 do begin
      bigdata(*,k2(i):k1(i)) = 0
   endfor

   ; Create parallel array of flags indicating where the good data are
   bigflag = bigdata
   bigflag(where(bigdata ne 0)) = 1

   k1 = (n*nobs+nacq+ntol)*npersec + koff  ; Array of start indexes of good data
   k2 = ((n+1)*nobs -ntol)*npersec + koff  ; Array of end indexes of good data

   ; Determine number of independent measurements within each block,
   ; and keep the minimum.
   nmeas = nobs - nacq - 2*ntol    ; Nominal number of measurements
   for i = 0, nblk-1 do begin
      istart = floor((k1(i)-koff)/float(lseq))
      iend   =  ceil((k2(i)-koff)/float(lseq))
      nmeas = nmeas < (iend - istart)
   endfor
   ; Create array to hold the averaged data
   avdat = fltarr(nant,lseq,nblk,nmeas)

   ; Loop over the 20 blocks (pointings) of data
   for i = 0, nblk-1 do begin

      ; Zero the accumulation arrays
      sum = fltarr(nant,lseq)
      nsum = intarr(nant,lseq)

      ; Determine start and end cycle numbers for this block.  In general
      ; both the start and end cycles will be partially full, but where
      ; the data are outside of the good data interval they have already
      ; been set to zero above, so they will not contribute to the sums.
      istart = floor((k1(i)-koff)/float(lseq))
      iend   =  ceil((k2(i)-koff)/float(lseq))<(istart+nmeas)

      ; Accumulate data for this block
      for j = istart,iend-1 do begin

         ; Determine start and end indexes of cycle
         j1 = j*lseq + koff
         j2 = j1 + lseq - 1

         ; Accumulate data, totalling data over the different cycles
         avdat(*,*,i,j-istart) = bigdata(*,j1:j2)
;         sum  =  sum + bigdata(*,j1:j2)
;         nsum = nsum + bigflag(*,j1:j2)
      endfor

      ; Determine the average for this block and store it
;      avdat(*,*,i) = sum/nsum

   endfor

   ; First determine good entries in observing sequence
   gdseq = where((*obseq.pharm)(0:lseq-1) ne -1)
   ; ...and good blocks in trajectory sequence (blocks 0,5,10, and 15
   ; are eliminated since antennas are slewing during these times.
   gdblk = [1,2,3,4,6,7,8,9,11,12,13,14,16,17,18,19]

   ; Select only the good values
   avdat = avdat(*,gdseq,gdblk,*)

   ; Perform fit to data contained in AVDAT, returning the peak level
   ; PK(nant,nf), the pointing offsets in degrees PO(nant,nf,2),
   ; the half-power-beamwidth HPBW(nant,nf,2), and rms error of the
   ; fit RMSERR(nant,nf,2)
   limfit,avdat,obseq,traj,pk,po,rmserr,debug=debug

   ; Restrict to 27-m antennas only
   nant = 2
   ; Convert PO from centidegrees to arcminutes
   po = po*60./100.

   pavg = (psig = fltarr(nant,2))
   ; Do a weighted average over frequency.  The weighting is
   ; just the Chi Square
   ;(a combination
   ; of the relative RMS error at each frequency and the peak flux value.)
   for i = 0, nant-1 do begin
      for j = 0, 1 do begin
         ;wt = reform((pk(i,*,j)/100.)*min(rmserr(i,*,j))/rmserr(i,*,j))
         wt = 1./rmserr(i,*,j)
         poff = reform(po(i,*,j))
         pavg(i,j) = total(poff*wt)/total(wt)
         psig(i,j) = sqrt(total((poff - pavg(i,j))^2)/n_elements(poff))
      endfor
   endfor

   fcfile = !defaults.workdir+'LIMPARM.DAT'
   ; Write out new HODO values into a file LIMPARM.DAT (overwrites any
   ; previously existing LIMPARM.DAT file).
   openw,fclun,/get_lun,fcfile
   printf,fclun,'/ LIMPOINT for ',traj.srcname+" "+tls.yrday+tls.timstr
   printf,fclun,pavg(0,0),psig(0,0),pavg(0,1),psig(0,1),$
      format='(F5.1,"(",F3.1,") HO,",F5.1,"(",F3.1,") DO")'
   printf,fclun,pavg(1,0),psig(1,0),pavg(1,1),psig(1,1),$
      format='(F5.1,"(",F3.1,") HO,",F5.1,"(",F3.1,") DO")'
   print,'/ LIMPOINT for ',traj.srcname+" "+tls.yrday+tls.timstr
   print,pavg(0,0),psig(0,0),pavg(0,1),psig(0,1),$
      format='(F5.1,"(",F3.1,") HO,",F5.1,"(",F3.1,") DO")'
   print,pavg(1,0),psig(1,0),pavg(1,1),psig(1,1),$
      format='(F5.1,"(",F3.1,") HO,",F5.1,"(",F3.1,") DO")'
   free_lun,fclun


   free_lun,lun

return
END