;+
; NAME:
;     APCALCHEK
; PURPOSE:
;     Main routine to analyze phase calibration data.  Responds to scans of type
;     PCAL (scan code 5).
; CATEGORY:
;     OVRO APC DATA ANALYSIS
; CALLING SEQUENCE:
;     apcalchek[,filename][,hrec,erec | ,after=after][,debug][,/cmdfile]
; INPUTS:
;     filename   the name of the file containing PCAL data.  If
;                  omitted, the file DAILY.ARC in directory
;                  !DEFAULTS.WORKDIR is assumed.
;     hrec       the record number of the PCAL 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:
;     erec       the end record number.  This is only needed as a
;                  convenience for the PROGRESS Meter, and may be omitted.
;     after      an optional time string of the standard form
;                  [yyyy.ddd ]hh:mm[:ss]m after which to start
;                  looking for a valid PCAL 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>\PCAL.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, get_centim, pldinit,
;     gaincor, pldout
; OUTPUTS:
; COMMENTS:
;     There will have to be modifications to this routine when
;     additional antennas are available.
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 22-Aug-1998 by Dale Gary (simple adaptation from PCALCHEK)
;     23-Aug-1998  DG
;       Fixed bug in polarization code determination.
;     25-Aug-1998  DG
;       Added check for new TRAJECTORY segments (but contents ignored
;       for now.)
;     05-Feb-2000  DG
;       Removed now unnecessary sin and cos lookup tables in GAINCOR call.
;-
pro apcalchek,filename,hrec,erec,after=after,debug=debug,cmdfile=cmdfile

   ; If the /cmdfile switch is set, open a "new" PCAL.MSG before doing
   ; anything.  If it exists it will be cleared.
   msglun = 0   ; Indicates PCAL.MSG does not exist
   IF (keyword_set(cmdfile)) THEN $
                  openw,msglun,/get_lun,!defaults.workdir+'PCAL.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

   ; EREC is just to let the progmeter know how many records there are.
   ; If not specified, use the value 100.
   IF (n_elements(erec) eq 0) THEN $
                  erec = hrec + 100

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

   ; Ask user which source was used.
   out = wstr('Phase Cal Source',text='Enter name of Phase Cal source',$
                              ['SOURCE NAME'],[''],['SRCNAME'])
   ; This is highly temporary--hard-code the source coordinates since
   ; source coordinates are not entered into the array parm segment yet.
   doy = lindgen(3) + tls.day
   CASE out.srcname OF
      '3C84': BEGIN
           ra = [49.9274,49.9274,49.9274]*240.D*1000.
           dec = [41.5025,41.5025,41.5025]*3600.D*1000.
         END
      '3C147': BEGIN
           ra = [85.6117,85.6117,85.6117]*240.D*1000.
           dec = [49.8474,49.8474,49.8474]*3600.D*1000.
         END
      '3C273': BEGIN
           ra = [187.2554,187.2554,187.2554]*240.D*1000.
           dec = [2.0627,2.0627,2.0627]*3600.D*1000.
         END
      ELSE: BEGIN
            if (result eq -1) then begin
               errstr = 'APCALCHEK: Unknown calibrator source specified.'
               IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                          ELSE ans = widget_message(errstr,/ERROR)
               return
            endif
         END
   ENDCASE

   ; Decode the header segment
   header = decode(data,!SEGM.HEADER)

   nant = header.nant

   ; Decode the config segment
   rec = rec + 1
   data = getdata(rec,a)
   cfg = decode(data,!SEGM.CONFIG)

   bx = cfg.bx(0:nant-1,0:nant-1)
   by = cfg.by(0:nant-1,0:nant-1)
   bz = cfg.bz(0:nant-1,0:nant-1)
   htdiff = cfg.htdiff(0:nant-1,0:nant-1)

   t = (doy - doy(0))*24.D*3.6D6

   ; Quadratic fit to coordinates -- returns polynomial coefficients
   ; A of form y(t) = A(0) + A(1)*t + A(2)*t^2.  [The REFORM()
   ; function is needed because otherwise poly_fit returns Array(1,3)
   ; instead of Array(3)--REFORM() converts Array(1,3) to Array(3)]

   a_ra  = reform(poly_fit(t,ra, 2))
   a_dec = reform(poly_fit(t,dec,2))

   ; Encapsulate it into a structure for passing to BDOTS and DELAY routines.
   headin = {bx:bx, by:by, bz:bz, htdiff:htdiff, a_ra:a_ra, a_dec:a_dec, year:tls.year, doy:doy}

   ; Decode the observing sequence
   rec = rec + 1
   data = getdata(rec,a)
   obseq = decode(data,!SEGM.OBSEQ)

   ; Check for trajectory segment
   rec = rec + 1
   data = getdata(rec,a)
   traj = decode(data,!SEGM.TRAJ)
   sz = size(traj)
   ; If the type of TRAJ is not 8 (not a structure), then TRAJ is not
   ; a trajectory structure, so backup on the record number.
   if (sz(n_elements(sz)-2) ne 8) then rec = rec-1

   ; Read in the channel offsets
   files = findfile(!defaults.dbdir+'*.off',count=n)
   offsets = replicate(0.0,25)
   if (n eq 0) then begin
      errstr = 'APCALCHEK: No channel offsets found.  Will use 0 for all.'
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
   endif else begin
      ; At least one file was found, so sort in numerical order and
      ; take the latest one.  Note that in general this is not the
      ; right thing to do--the closest one to the data's date should
      ; be used.
      files = files(sort(files))
      filename = files(n-1)
      openr,offlun,/get_lun,filename
      ; Read and discard the header lines
      hedlines = strarr(3)
      readf,offlun,hedlines
      ; Read the next line, and use it--really should use whichever is closest
      ; to data's date
      line = ''
      readf,offlun,line
      reads,line,offsets,format='(17X,25F5.1)'
      free_lun,offlun
   endelse

   ; The PCAL scan can have a variable number of entries.  Each cycle has
   ; a total of nsamp samples, each with nchan (=header.nws-nndw) channels,
   ; so declare an appropriately sized array.

   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 = obseq.nx                         ; Number of samples in one cycle
   tsample = indgen(nsamp)*cfg.sampintms    ; Time offsets [msec] corresponding
                                            ;   to each sample in one cycle
   bigdata = fltarr(nws,nsamp)              ; Big enough to hold one cycle
   nent = nws*nsamp                 ; Total number of entries in one cycle

   ; Obtain the polarization distribution
   ; NB: This code should be changed to check that ant 1 in at index 0 and
   ;      ant 2 at index 1.
   apol = bytarr(nsamp,8)    ; Create a poln list for each antenna
   for i = 0, nsamp-1 do apol(i,*) = cfg.antpol(0,*)    ; Set all polarizations to R
   apol(*,0) = cfg.antpol(*obseq.ppol and 1,0)     ; Replace ant 1 pol with obseq value
   apol(*,1) = cfg.antpol(ishft(*obseq.ppol and 4,-2),1)  ; Replace ant 2 pol with obseq value
   pol = bytarr(nws,nsamp)     ; Parallel array to bigdata, containing poln code
   for j = 0, nant-1 do begin
      for i = 0, j do begin
         k = header.lchoff(i,j)
         if (i eq j) then begin
            pol(k,*) = apol(*,i)
         endif else begin
            ; Channel poln will always be lesser of two (only possibilities
            ; are 0-0, 0-3, 1-1, 1-3, or 3-3, corresp. to polns 0 0 1 1 3)
            pol(k,*) = apol(*,i) < apol(*,j)
            pol(k+1,*) = pol(k,*)
         endelse
      endfor
   endfor

   ; List of indexes where noise diode is off (OFFIDX)
   ; as obtained from observing sequence (OBSEQ structure)
   offidx = where(*obseq.pnd EQ 2)
   nfp = n_elements(offidx)   ; Number of frequencies/polarizations

   recn = 1                         ; Record number of current data record
   j = 0
   rec = rec + 1
   ndat = header.nws-nndw           ; Number of words of data in record

   ; Capture end of file error to indicate that we reached the end of the
   ; file.
   ON_IOERROR,done
   nnow = 0
   nval = 0

   ; Turn off math exception handling to eliminate annoying report of
   ; Floating illegal operand every time a NaN is processed.
   !EXCEPT = 0

   ; Initialize PLD file and get its logical unit number (pldlun)
   result = pldinit(pldlun,'APCALDAT',header,cfg,obseq)

   ; Start the progress meter
   id = progmeter(/init,buttontext="Abort")
   IF (result ne 0) THEN BEGIN
      errstr = 'APCALCHEK: PLD file could not be created.'
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
      return
   ENDIF

   ; Number of samples before start of DOSEQ
   koff = 2

   ; Read records until an error occurs or some other condition causes
   ; us to stop.
   while (1) do begin

      datoff = header.rdoff+koff*nws  ; Location of first sample in data
                                      ; koff to skip first two samples
      ; Read data from a record
      data = getdata(rec,a)

      ; Update progress meter.  Note that for a growing file, erec cannot
      ; be known so if the record number exceeds the estimated erec then
      ; double the estimate.
      if (rec gt erec) then erec = hrec + (erec-hrec)*2
      if(progmeter(id,float(rec-hrec)/(erec-hrec)) eq 'Cancel') then goto,done

      if(data(2) eq 0) then begin

         ; If data record is all zeros, we must have gone beyond the valid
         ; data, so wait for 1 s and try again.
         wait,1

      endif else begin

         ; This is good data, so decode the time/label field
         tls = tl_decode(data)
         if (tls.segmentcode ne !SEGM.DATA) then goto,done
         if (nval eq 0) then begin
            tcycle = tls.msec
         endif
         nval = nws*(nsr-koff)             ; Number of data values in a record

         koff = 0    ; Reset koff after first record has been read

         ; Determine how many more data values we need from DATA array to
         ; fill BIGDATA array.

         inner:
           n = nent - nnow     ; Number of entries yet to read
           nmore = nval < n    ; Smaller of nval and n (nval is number in one record)
           bigdata(nnow) = data(datoff:datoff+nmore-1) ; Convert 1d DATA array to 2d BIGDATA array
           nnow = nnow + nmore ; Advance pointer to current entry
           if (nnow eq nent) then begin

              ; BIGDATA is full, so time to prepare it for writing to disk

              flagool,bigdata,header,obseq   ; Flag data for out of lock

              ; Correct for complex gain and atten
              smldata = gaincor(bigdata,pol,header,cfg,obseq,tcycle,headin,offsets)

              ; GAINCOR returns an array with all 3 poln states, but since APCALCHEK only
              ; takes I data, we get rid of the other two.
              smldata = reform(smldata(*,*,2))
              result = pldout(pldlun,tcycle,smldata)    ; Write data to plot-data-file

              nnow = 0
              datoff = datoff + nmore         ; Offset to next data in record
              nval = nval - nmore             ; Reset number of values to go in record
              tls = tl_decode(data)
              tcycle = tls.msec+(nmore/header.nws)*cfg.sampintms
              if (nval ne 0) then begin
                 goto, inner
              endif
              ; nval is zero, so we need to get the time from the next record
           endif

         rec = rec + 1
      endelse
   endwhile

done:
   ON_IOERROR,NULL
   ; Turn on math exception handling again
   !EXCEPT = 1
   ; Destroy progress meter
   status = progmeter(id,/destroy)
   result = pldout(pldlun,-2L,smldata*0-2.0)
   free_lun,pldlun
   free_lun,lun

return
END