;+
; NAME:
;     GCALCHEK
; PURPOSE:
;     Main routine to analyze gain calibration data.  Responds to
;     scans of type GCAL (scan code 4) or TEST DATA (scan code 1).
; CATEGORY:
;     OVRO APC DATA ANALYSIS
; CALLING SEQUENCE:
;     gcalchek[,filename][,hrec | ,after=after][,gains][,debug][,/cmdfile]
; INPUTS:
;     filename   the name of the file containing GCAL data.  If
;                  omitted, the file DAILY.ARC in directory
;                  !DEFAULTS.WORKDIR is assumed.
;     hrec       the record number of the GCAL 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]m after which to start
;                  looking for a valid GCAL 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>\GCAL.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, lasthrec, newscan, getdata, tl_decode, decode, gcfit, nlfit,
;     write_off, write_gcparm, write_gainrec
; OUTPUTS:
;     gains      optional output parameter to return the receiver, noise
;                  diode, and rms spectra in an array of size GAINS(3,NANT,86)
;                  or GAINS(3,NANT,172), depending on the epoch
; COMMENTS:
;     Starting on 2003.310 (5-Nov-2003), the GAINCAL consists of both
;     RCP and LCP measurements.  This routine should work for epochs both
;     before and after this.
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 17-Jan-1998 by Dale Gary
;     08-Mar-1998  DG
;       Updated to implement a somewhat different scheme.  The
;       nonlinearity parameter is determined only from 0 and 5 DB
;       data, in routine NLFIT.
;     15-Mar-1998  DG
;       Added cmdfile switch for unattended operation.
;     18-Jun-1998  DG
;       Added koff = 2 as number of samples before DOSEQ starts.
;       Without this, we were off by one harmonic.
;     21-Jun-1998  DG
;       Changed from reading a(rec) to getdata(rec,a) throughout,
;       so that routine should work under UNIX.
;     25-Aug-1998  DG
;       Added check for new TRAJECTORY segments (but contents ignored
;       for now.)
;     23-Feb-1999  DG
;       Convert to use NEWSCAN.
;     07-Apr-1999  DG
;       Implement early form of GCAL.DIAG file, to write out diagnostics
;       for use by DIAGNOSE routine.
;     01-Dec-1999  DG
;       Changed to write the GPARM.REC file (and write the gain parameter
;       segment directly into the .ARC file if appropriate).  This reflects
;       changes in WRITE_GAINREC.PRO.
;     11-Jan-2000  DG
;       Eliminate hardwired directory
;     16-Dec-2000  DG
;       Added optional argument, GAINS that, if present, causes GCALCHEK
;       to skip all output and return the receiver and noise diode values
;       to the calling program instead.  This will be used to go through
;       a large number of files and accumulate receiver and noise diode
;       statistics.
;     11-Mar-2001  DG
;       Made sin/cos channel offsets consistent by enforcing that all AB
;       channels have same offset, all BC channels, etc.
;     06-Nov-2003  DG
;       Several changes to make it work for new RCP/LCP measurements.  The
;       gain parameters are calculated only for RCP (as before), but the
;       GAINS array now contains both polarizations.
;     09-Jan-2004  DG
;       Slight changes to plot_gcparm and gcfit calls to use polarization
;       keyword.
;     18-Jul-2004  DG
;       Significant changes to use NAN as non-data values instead of 0 or -1.
;     23-Aug-2005  DG
;       Fix bug handling older GCAL data (due to GCSTAT2 undefined).
;-

pro gcalchek,filename,hrec,gains,after=after,debug=debug,cmdfile=cmdfile

   common GC_prgmtr, id

   id = progmeter(/init,label='Analyzing Gain Calibration')
   gains = 0

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

   res = progmeter(id,0.1)

   ; Read header info, etc.
   rec = newscan(a,nrec,hrec,header,obseq,cfg,traj,refcal,gparm)

   if (rec eq -1) then return  ; Bad scan

   res = progmeter(id,0.2)

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

   nant = header.nant
   tls = header.tls

   ; The GCAL should have obseq.nx*obseq.nt (178*28) entries, each with
   ; nchan (=header.nws-nndw=25) channels, so declare an appropriately
   ; sized array.

   nndw = header.ndbw + header.nfgw + 1     ; Number of non-data words
   nchan = header.nws-nndw                  ; Number of channels
   nbl  = (nchan-nant)/2                    ; Number of baselines
   nsamp = obseq.nx*obseq.nt                ; Number of samples (178*28)

   ; Get index of active antennas IANTS (e.g. [1,2] => [0,1], [4,5,6] => [2,3,4]
   antlist = [1,2,4,5,6,7,8]   ; Translation from antenna number to index
   aa = header.aatab(where(header.aatab ne 255))
   iants = intarr(nant)
   for i = 0, nant-1 do begin
      iants[i] = where(antlist eq aa[i])
   endfor
   ; Nant x nant matrix of channel offsets and MUXARR.  Also FG pointer array.
   muxarr = (choff = intarr(nant,nant))
   for i = 0, nant-1 do begin
      for j = 0, nant-1 do begin
         muxarr[i,j] = cfg.muxarr(iants[i],iants[j])
         choff[i,j]  = header.lchoff(iants[i],iants[j]) - nndw
      endfor
   endfor
   ; Form array representing the indexes of upper part of matrix
   upper = intarr(nbl)
   k = 0
   for i = 0, nant-2 do begin
      for j = i+1, nant-1 do begin
         upper[k] = i*nant+j
         k = k + 1
      endfor
   endfor

   bigdata = intarr(nchan,nsamp) + '7FFF'x  ; Big enough to hold all data (init to SATURATED)

   ; The number of records to expect should be nrecs=nsamp/header.nsr
   ; (178*28/33 = 151).  Note that there should be one sample in the 152nd
   ; record, but CPC doesn't send it.
   nrecs = nsamp/header.nsr
;   rec = rec + 1
   datoff = header.rdoff+nndw
   ndat = header.nws-nndw

   ; Read data from the file into a single large array, BIGDATA, of
   ; size (25,178*4*7), which is just a long list of all 25 data channels
   ; Loop over the 151 records
   for i = 0,nrecs-1 do begin
      data = getdata(rec,a)                      ; Read a record
      tls = tl_decode(data)
      case tls.segmentcode of
         !SEGM.DATA: begin
               ; 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)
               endfor
             end
         !SEGM.EOS: goto,bail
         !SEGM.HEADER: goto,bail
         ELSE:
      endcase
      rec = rec + 1                      ; Increment record number
   endfor

bail:
   res = progmeter(id,0.3)
   free_lun,lun

   nasteps = 7 ; Number of attenuation steps (0,5,10,15,20,35,20*DB)
   nrpt = 4    ; Number of times entire sequence should repeat
   IF (nasteps*nrpt NE obseq.nt) THEN BEGIN
      errstr = ['GCALCHEK: Problem with observing sequence.',$
                  'Number of cycle repeats is different than expected.']
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
      return
   ENDIF

   ; List of indexes where noise diode is off (OFFIDX) and on (ONIDX)
   ; as obtained from observing sequence (OBSEQ structure)
   offidx = where(*obseq.pnd EQ 2)
   onidx  = where(*obseq.pnd EQ 1)
   IF (n_elements(offidx) NE n_elements(onidx)) THEN BEGIN
      errstr = ['GCALCHEK: Problem with observing sequence.',$
                           'Noise diode duty cycle not 1:1.']
      IF (keyword_set(cmdfile)) THEN printf,msglun,errstr,format='(a)' $
                                ELSE ans = widget_message(errstr,/ERROR)
      return
   ENDIF

   ; Declare storage for main GCAL array, GCTEMP, a 5-dimensional array of
   ; size (2,7,25,86,4) where:
   ;    first  index     0..1  ND off/on
   ;    second index     0..6  Attenuator state
   ;    third  index     0..24 Channel number
   ;    fourth index     0..85 Frequency
   ;    fifth  index     0..3  Measurement interation
   gctemp = fltarr(2,nasteps,nchan,n_elements(offidx),nrpt)

   koff = 2   ; Number of samples before DOSEQ starts
   ; Apportion data in BIGDATA array into GCTEMP. Note that data that
   ; are expected to be bad are flagged with zeroes in the *OBSEQ.PND array.
   ; Thus, these entries do not appear in either the OFFIDX or ONIDX list
   ; of indexes, and therefore are excluded from the GCTEMP array.
   for k = 0, nrpt-1 do begin                ; Loop over the 4 measurements
      for i = 0, nasteps-1 do begin          ; Loop over the 7 atten states
         ioff = offidx+i*obseq.nx          ; Shift off and on indexes by
         ion  =  onidx+i*obseq.nx          ;   length of observing sequence
         gctemp(0,i,*,*,k) = bigdata(*,koff+ioff+k*nasteps*obseq.nx)  ;ND off state
         gctemp(1,i,*,*,k) = bigdata(*,koff+ion +k*nasteps*obseq.nx)  ;ND on  state
      endfor
   endfor

   res = progmeter(id,0.4)

   ; Save a copy of GCTEMP and set all saturated data ('7FFF'X) to NAN
   ; so that it does not affect total.
   gcgood = gctemp
   saturated = where(gctemp eq '7FFF'X,nsat)
   IF (nsat NE 0) THEN gcgood[saturated] = !values.f_nan

   ; Calculate the sum over the 5th index (i.e. the 4 measurements) and
   ; the sum of the squares GCTOT2.  The GCTOT and GCTOT2 arrays have
   ; size (2,7,NBL,NFP) where NBL = number of baselines, NFP = number of
   ; frequencies*poln
   gctot = total(gcgood,5,/NAN)
   gcsqr = gcgood^2.
   gctot2 = total(gcsqr,5,/NAN)

   ; GCFLAG is a parallel flag array to GCGOOD.  Where GCGOOD is good,
   ; FLAG=1, and where GCGOOD is saturated, GCFLAG=0.  This is handy
   ; for indicating good vs bad data, but also for summing total numbers
   ; of good data.  This general technque is used often below and in the
   ; GCFIT routine.
   gcflag = fix(gcgood*0) + 1
   IF (nsat NE 0) THEN gcflag[saturated] = 0

   ; Use the GCFLAG array to determine the total number of good values, so
   ; that it can be used for normalization.  Divide the total by the sum
   ; to get the average (when gcsum is zero (no good data) we set it to one
   ; since gcavg will be zero anyway by virtue of GCTOT being zero there).
   gcsum = fix(total(gcflag,5))
   gcavg = gctot/(gcsum>1)

   ; Calculate the variance
   gc2 = gcgood
   for i = 0, nrpt-1 do begin
      gc2[*,*,*,*,i] = (gcgood[*,*,*,*,i] - gcavg)^2.
   endfor
   gcvar = total(gc2,5,/NAN)/(gcsum>1)

   ; Save the average values back into GCTEMP, then flag GCTEMP values to NAN
   ; wherever GCVAR exceeds 10^2. GCTEMP now has size (2,7,25,86)
   gctemp = gcavg
   bad = where(gcvar GT 10^2 OR gcsum EQ 0,nbad)
   IF (nbad NE 0) THEN gctemp[bad] = !values.f_nan

   ; Also flag values for which GCSUM equals zero, if any
;   bad = where(gcsum EQ 0)
;   IF (bad(0) NE -1) THEN gctemp(bad) = -1

   ; Calculate the 25 channel offsets by using the 6th attenuation
   ; setting (35 DB).  GCOFF and the arrays derived from it now have
   ; dimensions (2,25,86)
   gcoff = reform(gctemp[*,5,*,*])

   ; Convert to a flag array to use for normalization, as before.
   gcofflg = fix(gcoff*0) + 1
   bad  = where(finite(gcoff) EQ 0,nbad)
   IF (nbad  NE 0) THEN gcofflg[bad] = 0

   if (n_params() lt 3) then begin
      ; Diagnostics: Open DIAG file
      OPENW,dlun,/get_lun,!DEFAULTS.DBDIR+'GCAL.DIAG'
      ; Write date/time string as first line of DIAG file
      PRINTF,dlun,'GCAL '+tls.yrday+strmid(tls.timstr,0,9)+' UT'
   endif

   ; Finally, determine the offsets for each of the 25 channels.
   ; The inner TOTAL function sums over ND on/off, the outer one
   ; sums over frequency.
   offsets = total(total(gcoff,1,/NAN),2,/NAN)/((total(total(gcofflg,1),2))>1)
   oldoff = offsets[0:nant-1]
   if (keyword_set(debug)) then $
         print,'Channel Offsets: '+string(oldoff,format='(8F6.1)')

   ; Reduce consideration to only the 5 total power channels, and
   ; subtract offsets.  Bad data remain flagged with NAN
   gcdata = gctemp[*,*,0:nant-1,*]
   gcflag = fix(gcdata*0) + 1
   bad = where(finite(gcdata) eq 0,nbad)
   if (nbad ne 0) then gcflag[bad] = 0

   ; Now ready for loop over different nonlinearity params
   abest = replicate(0.0,nant)
   resids = fltarr(nant,101)

   res = progmeter(id,0.5)

   ; Calculate offsets using zero nonlinearity
   ; As of 2003.310, the Gaincal does both R and L polarization, but we
   ; only need offsets for one polarization, so just use the first 86 entries
   gcparm = gcfit(gcdata[*,*,*,0:85],gcflag[*,*,*,0:85],abest,resid,oldoff,/refine_off,msglun=msglun,pol='R')
   newoff = oldoff

   ; Adjust channel offsets for consistency (all AB offsets should be the same, all
   ; BC offsets, etc.)

   all = choff[upper]
   for i = -1, 1 do begin
      ilist = all[where(muxarr[upper] eq i)]
      offsets[ilist  ] = avg(offsets[ilist  ])
      offsets[ilist+1] = avg(offsets[ilist+1])
   endfor

   ; Determine backend offsets from these psuedo-channel offsets
   ab = (all[where(muxarr[upper] eq 0)])[0]
   bc = (all[where(muxarr[upper] eq -1)])[0]
   ca = (all[where(muxarr[upper] eq 1)])[0]
   ab = [ab,ab+1]
   bc = [bc,bc+1]
   ca = [ca,ca+1]
   beoff = [offsets[0:5],offsets[ab],offsets[bc],offsets[ca]]

   ; As of 2003.310, the Gaincal does both R and L polarization, but we
   ; only calculate nonlinearity for one polarization, so just use the first 86 entries
   ; Note that the LNA is different for the two polarizations, and could contribute
   ; to the overall nonlinearity, but we ignore this
   ; Get the nonlinearity for the five antennas (an array of size 5)
   if(keyword_set(debug)) then nl = nlfit(gcdata[*,*,*,0:85],gcflag[*,*,*,0:85],newoff,resid,/debug) $
                          else nl = nlfit(gcdata[*,*,*,0:85],gcflag[*,*,*,0:85],newoff,resid)

   res = progmeter(id,0.6)

   nfp = n_elements(gcdata[0,0,0,*])
   gcstat = strarr(3,5,nant)
   s2nspec = fltarr(nant,nfp,2)
   if (n_params() lt 3) then begin
      antstr = ['1','2','4','5','6','7','8']
      typstr = ['IF','Noise Diode','Gain','Nonlinearity','Offsets']
      for i = 0, nant-1 do begin

         ; Loop over data type, to initialize status strings
         for j = 0, 4 do begin
            gcstat(0,j,i) = "Ant "+antstr[i]+" "+typstr[j]+" status:"
            gcstat(1,j,i) = ' '
            gcstat(2,j,i) = '    '
         endfor

         ; Diagnostic - S/N ratio from the 4 measurements
         ; Loop over ND off/on
         for j = 0, 1 do begin
            s2nspec(i,*,j) = gcavg(j,1,i,*)/sqrt(gcvar(j,1,i,*)>1)
;         plot,s2nspec(i,*,j),psym=-1
            s2n = total(s2nspec(i,*,j))/total(gcvar(j,1,i,*) gt 0)
            ; Show good S/N if > 100, marginal if > 50, low if < 50
            CASE 1 OF
              s2n gt 100.: gcstat(1,j,i) = 'Good S/N '+string(fix(s2n),format='(I4)')+':1'
              s2n gt  50.: gcstat(1,j,i) = 'Marg S/N '+string(fix(s2n),format='(I4)')+':1'
              s2n le  50.: gcstat(1,j,i) = 'Low  S/N '+string(fix(s2n),format='(I4)')+':1'
            ENDCASE
         endfor

         if (abs(nl(i)) ge 0.6) then gcstat(2,3,i) = 'WARN: Very Nonlinear'
         if (abs(nl(i)) ge 1.0) then gcstat(2,3,i) = 'BAD : Max Nonlinearity'
      endfor
   endif

   res = progmeter(id,0.65)

   ; Calculate the attenuator values, and determine the receiver
   ; and noise diode increment spectra.  This call will create a PLD
   ; file containing the receiver level, noise diode increment, and
   ; residual spectra.

   if (n_params() eq 3) then begin
      ; Just determine the IF levels and return to caller
      gcparm = gcfit(gcdata[*,*,*,0:85],gcflag[*,*,*,0:85],nl,resid,newoff,tls,gcstat,s2nspec[*,0:85,*],msglun=msglun,gains,/ifonly,pol='R')
      if (tls.yrday ge 2003.310) then begin
         gcparm2 = gcfit(gcdata[*,*,*,86:*],gcflag[*,*,*,86:*],nl2,resid2,newoff2,tls,gcstat2,s2nspec[*,86:*,*],msglun=msglun,gains2,/ifonly,pol='L')
         gains = [[[gains]],[[gains2]]]
         dims = size(gains,/dim)
         gains = reform(gains,dims[0],dims[1],dims[2]/2,2)
      endif
      res = progmeter(id,/destroy)
      return
   endif else begin
      gcparm = gcfit(gcdata[*,*,*,0:85],gcflag[*,*,*,0:85],nl,resid,newoff,tls,gcstat,s2nspec[*,0:85,*],msglun=msglun,gains,pol='R')
      if (tls.yrday ge 2003.310) then begin
         nl2 = nl
         newoff2 = newoff
         gcstat2 = gcstat
         gcparm2 = gcfit(gcdata[*,*,*,86:*],gcflag[*,*,*,86:*],nl2,resid2,newoff2,tls,gcstat2,s2nspec[*,86:*,*],msglun=msglun,gains2,pol='L')
         gains = [[[gains]],[[gains2]]]
         dims = size(gains,/dim)
         gains = reform(gains,dims[0],dims[1],dims[2]/2,2)
      endif
   endelse

   res = progmeter(id,0.7)

   ; Enter corrected total power offsets back into offsets array
   offsets(0:nant-1) = newoff

   if (keyword_set(debug)) then begin
      print,'Revised Offsets: '+string(newoff,format='(8F6.1)')
      print,' 5 DB corrctn factor: ',reform(gcparm(2,*))
      print,'10 DB corrctn factor: ',reform(gcparm(3,*))
      print,'20 DB corrctn factor: ',reform(gcparm(4,*))
      print,'Nonlinearity factors: ',nl
   endif

   ; Write information to standard database files.  These routines
   ; check whether data for this date and time already exist in the
   ; files, and if so will ask for confirmation to overwrite.
   result = write_off(offsets,tls,msglun=msglun)
   result = write_gcparm(gcparm,tls,msglun=msglun)
   noupdate = write_gainrec(header,gcparm,offsets,gains,msglun=msglun,outrec)

   ; Determine whether OUTREC has 3 or 6 records
   n_rec = n_elements(outrec[0,*])

   ; Close the file and reopen it for writing
   openu,lun,/get_lun,filename
   a = assoc(lun,intarr(1024))
   ; Get the record after the header record of the GCAL scan, and see if it
   ; has a segmentcode of zero or !SEGM.GPARM (indicates that it can be written to)
   rec = hrec
   ; Loop over the three records of the GAIN PARAMETER segment
   for i = 0,n_rec-1 do begin
      ; Get a record from the .ARC file
      data = getdata(rec,a)
      tlsnew = tl_decode(data)
      ; Check that its segmentcode is as expected
      if (tlsnew.segmentcode eq 0 or tlsnew.segmentcode eq !SEGM.GPARM) then begin
         outrec[0,i] = tlsnew.nrc  ; Put record number into OUTREC (from WRITE_GAINREC)
         a(rec) = outrec[*,i]      ; Write OUTREC to the file in place of the current one
         rec = rec + 1             ; Increment to next record
      endif else begin
         ans = dialog_message(['GCALCHEK: GCAL does not have room for ',$
                               'GAIN PARAMETER segment.'],/error)
         goto, skip     ; Wrong segmentcode so bail
      endelse
   endfor

   ; Write the data out to the .REC file unless NOUPDATE = 1 (NOUPDATE is returned
   ; by WRITE_GAINREC, and indicates whether this is a current GCAL or an older
   ; one determined off-line).
   if (noupdate eq 0) then begin
      ;N.B. The record number will be used by ARCHIVE to as the rec #
      ; of GAINPARM segment to write into the index record, for subsequent
      ; headers received by ARCHIVE
      ON_IOERROR, skip             ; If the directory doesn't exist, skip writing the file
      openu,reclun,/get_lun,!defaults.ephemdir+'gparm.rec'
      b = assoc(reclun,intarr(1024))
      for i = 0, n_rec-1 do b[i] = outrec[*,i]
      free_lun,reclun
   endif

skip:
   free_lun,lun
   IF (keyword_set(cmdfile)) THEN free_lun,msglun

   res = progmeter(id,0.8)

   plot_gcparm,tls,nant,gcstat,pol='r'
   if (n_elements(gcstat2) gt 0) then plot_gcparm,tls,nant,gcstat2,pol='l'

   res = progmeter(id,/destroy)

   ; Write out Diagnostics info and close file
   if (n_elements(gcstat2) gt 0) then PRINTF,dlun,gcstat2,format='(a)' else PRINTF,dlun,gcstat,format='(a)'
   FREE_LUN,dlun

return
END