;+
; NAME:
;     DAILY
; PURPOSE:
;     Automatic analysis and solution of PCAL data in a given file.
; CATEGORY:
;     OVSA APC CALIBRATION ANALYSIS
; CALLING SEQUENCE:
;     refscan = daily([filename][,ref = refscan | ,iref=iref][,/update] $
;                     [,/noprinter][,/nowrite],[/interactive][,sn=sn])
; INPUTS:
;     filename    The name of the .ARC file containing the data to be
;                   analyzed.  If omitted, the user is prompted for a filename.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     ref         A scan structure as returned from a previous call to DAILY.
;                   This is used when two different files are to be fit using
;                   a single reference scan.  Overrides keyword IREF.
;     iref        Specifies which scan is to be used a the reference--overrides
;                   the automatic selection that is done based on S/N.  This
;                   keyword is ignored if REF keyword is used.
;     update      A switch that, if set, causes scans that have already been
;                   analyzed not to be reanalyzed.
;     noprinter   By default, an overview page is printed on a local printer.
;                   This keyword causes the results to be shown on the display
;                   instead.
;     nowrite     By default, the results are written to the data (.ARC) file.
;                   This keyword causes the writing to be skipped.
;     interactive By default, the analysis is done automatically.  This
;                   keyword allows the intermediate results of the fitting
;                   to be adjusted interactively.
;     sn          Specifies Signal to Noise ratio to use in WPHZFIT/PHZFIT.  By
;                   default, this is 2.0, but can be lowered to 1.5 or so if too
;                   many baselines are missing due to noisy data.
; ROUTINES CALLED:
;     openarc, getdata, find_index, decode, analyze, amphit, cal_info, plot_phz,
;     daily_plot
; OUTPUTS:
;     refscan    A scan structure containing the data used for the reference
;                  scan, for use in a subsequent call to DAILY.
; COMMENTS:
;     Generally the routine will be called once without any reference scan
;     specified.  Sometimes, however, the data for one day is split into two
;     separate files.  In that case, call DAILY on the first file, without
;     specifying a reference scan, then use the output reference scan as
;     input for a second call, to the other file.  (Another, better way to
;     do this is to combine the two files using CONCAT_OVSA_FILE.)
; SIDE EFFECTS:
;     A page is printed to the local printer, unless overridden by the
;     /NOPRINTER keyword.  The results are written to the data file unless
;     overridden by the /NOWRITE keyword.
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 21-May-2000 by Dale E. Gary
;     29-May-2000  DG
;       Basically in final form
;     01-Jun-2000  DG
;       Added REFSCAN output for handling two separate files with a single
;       reference scan.
;     17-Jul-2000  DG
;       Several changes to make it honor both slope and offset of fits,
;       and also fix a bug in which I had not been correcting for the
;       fits before averaging the daily cals.  Change output time to be
;       a reference time to center of scan.  Add UPDATE switch.
;     23-Jul-2000  DG
;       Added IREF keyword, to give some control over which scan is used
;       as the reference.  Added NOPRINTER, NOWRITE, and INTERACTIVE switches.
;       Implemented WRITE_DAILY routine to write results to the .ARC file.
;     27-Jul-2000  DG
;       Fixed a couple of small bugs introduced in above changes.
;     06-Sep-2005  DG
;       Added changes to work with 7 antennas (I thought I had already made
;       these changes, but I do not see them...)
;     06-Nov-2005  DG
;       Phase closure does not work unless phases are corrected for baseline
;       ORDER and backend phase BEPHZ.  I wrote a new routine CAL_INFO that
;       is used to correct the phases for use in WPHZFIT, where the phase
;       closure is done.  The reference phases also are corrected for
;       comparison with the DAILY phase results.  Finally, the values are
;       converted back for insertion into the data file with WRITE_DAILY()
;     18-Jan-2006  DG
;       Added SN keyword to allow lowering of signal to noise limit, in
;       case of really crappy data.  SN defaults to 2, but can be lowered
;       to 1.5 or so if there are too many missing baselines.
;     27-Jul-2010  DG
;       I found a couple of commented ilnes that were not supposed to be
;       commented--applying order and bephz to the reference calibration.
;       These are no uncommented.
;     15-May-2014  DG
;       Added code to properly skip a scan when "skip" is selected in wphzfit call.
;       The code was previously not doing the proper thing in this case.
;     16-May-2014  DG
;       Changed call to WPHZFIT to include source name in the window title.
;     30-May-2014  DG
;       Added some printing of information on scans to be processed, so that
;       the user has a better idea of why sometimes scans are skipped.
;-

function best_scan,scan

   ; This routine will examine S/N and existence of frequencies and select
   ; the best scan to use as a reference.  The output is simply the index of
   ; the best scan.
   nscan = n_elements(scan)
   junk = where(scan[0].amf[0,0,*] ne 0,nbl)
   sn = reform(scan.amf[1,*,0:nbl-1]/scan.amf[2,*,0:nbl-1])
   snmax = sn[*,*,0]
   iref = 0
   for i = 1, nscan-1 do begin
      if (total(sn[*,*,i] gt snmax) gt n_elements(snmax)/2.) then begin
         snmax = sn[*,*,i]
         iref = i
      endif
   endfor

return,iref
end

function daily, filename, ref=refscan, update=update, iref=iref, noprinter=noprinter, $
                nowrite=nowrite, interactive=interactive, sn=sn

   if (not keyword_set(sn)) then sn = 2.0

   ; If no filename was specified, ask the user for a file.  If the user cancels,
   ; exit.
   if (n_elements(filename) eq 0) then begin
      nofile = 1
ask:
      filename = dialog_pickfile(path=!defaults.datadir,filter='*.arc',$
                    title='Select file containing PCAL measurements')
      if (filename eq '') then return,0
   endif

   lun = openarc(filename,a,nrec)

   ; If no file could be opened, and a filename was specified in the
   ; command line, return quietly.  If no filename was specified in the
   ; command line, then we are in interactive mode and it is okay to
   ; ask for another file.
   if (lun eq 0) then if (n_elements(nofile) eq 0) then return,0 else goto,ask

   ; Check for existence of an INDEX record
   idx = ovsa_get_index(a,nrec)

   ; Create a flag to indicate presence of valid index structure
   if (idx.nscans eq 0) then noidx = 1 else noidx = 0

   if (noidx) then begin

      ; There is no INDEX, so we have to find the information the hard way...

      ; Read the ACAL headers from the file
      seglist = get_segments(a,nrec)
      nadded = findsegentries(seglist,ihed,iseg)
;      seglist = seglist[ihed]
      iacal = where(strpos(seglist,'Header: PHASE CAL') ne -1,nacal)
      if (nacal ne 0) then seglist = seglist[iacal]
   endif else begin
      iacal = where((*idx.pscan).scancode eq !scan.pcal,nacal)
   endelse

   if (nacal eq 0) then begin
      ; If a filename was specified, we are not in interactive mode, so
      ; exit quietly.
      if (n_elements(nofile) ne 0) then begin
         print,'DAILY: File contains no PCAL scans.'
         return,0
      endif

      ans = dialog_message(['DAILY: File contains no PCAL scans.','Please select another file.'],/error)
      free_lun,lun
      goto,ask
   endif

   nanal = 0
   nf = 6
   nblmax = 11

   ; Fill in scanlist structure
   scan = replicate({sel:' ',srcname:'',date:'', time:'',srec:0,erec:0,anal:' ',$
                     tls:get_tl_struct(),nant:0, amf:fltarr(6,nf,nblmax), poff:fltarr(nblmax), pslp:fltarr(nblmax)},nacal)
   if (noidx) then begin
      ; No INDEX, so go through the file one scan at a time

      ; Loop over all scans
      for i = 0, nacal-1 do begin
         scan[i].srec = fix(strmid(seglist[i],13,5))
         scan[i].erec = fix(strmid(seglist[i],19,5))

         scan[i].sel  = '*'
         ; Read the header
         hrec = scan[i].srec
         data = getdata(hrec-1,a)   ; Read the header
         header = decode(data,!SEGM.HEADER)
         scan[i].tls = header.tls
         scan[i].nant = header.nant
         scan[i].date = header.tls.date
         scan[i].time = strmid(header.tls.timstr,1,8)

         ; Find the trajectory segment record that corresponds to this scan
         rec = findseg(a,hrec,nrec,!SEGM.TRAJ,idx=idx)
         ; Decode the trajectory segment, if any
         data = getdata(rec-1,a)
         traj = decode(data,!SEGM.TRAJ)
         scan[i].srcname = traj.srcname
      endfor
   endif else begin
      ; There is an INDEX, so get the info from the INDEX structure
      scan.srec = (*idx.pscan)[iacal].srec
      scan.erec = (*idx.pscan)[iacal].erec
      ; Read the header
      for i = 0, nacal-1 do begin
         hrec = scan[i].srec
         data = getdata(hrec-1,a)   ; Read the header
         header = decode(data,!SEGM.HEADER)
         scan[i].nant = header.nant
         scan[i].tls  = header.tls
      endfor

      scan.sel  = '*'

      scan.srcname = (*idx.pscan)[iacal].srcname
      for i = 0, nacal-1 do begin
         datarr = cvdoy((*idx.pscan)[iacal[i]].syear,(*idx.pscan)[iacal[i]].sday)
         scan[i].date = string(datarr,format='(I4,"/",I2.2,"/",I2.2)')
         scan[i].time = strmid(msec2str((*idx.pscan)[iacal[i]].smsec),1,8)
      endfor
   endelse

   ; Now have the PCAL scan list, but some may not be the correct sequence (51 DOSEQ),
   ; so go through and check each one
   for i = 0, nacal-1 do begin
      hrec = scan[i].srec
      rec = findseg(a,hrec,nrec,!SEGM.OBSEQ,idx=idx)
      data = getdata(rec-1,a)
      obseq = decode(data,!SEGM.OBSEQ)
      if (obseq.id ne 51) then begin
         scan[i].sel = ''  ; Deselect if this scan doesn't use 51 DOSEQ
         print,scan[i].srcname, scan[i].date, scan[i].time, scan[i].srec, scan[i].erec, format='(A8," scan at ",A10,A8,1X,I5,":",I5," will be skipped (wrong obseq)")'
      endif else begin
         print,scan[i].srcname, scan[i].date, scan[i].time, scan[i].srec, scan[i].erec, format='(A8," scan at ",A10,A8,1X,I5,":",I5," will be used")'
      endelse
   endfor

   free_lun,lun

   ; Have to read the number of antennas from the first scan
   lun = openarc(filename,a,nrec)
   data = getdata(scan[0].srec-1,a)
   header = decode(data,!SEGM.HEADER)
   ants = where(header.aatab gt 2 and header.aatab lt 255,nant)
   free_lun,lun

   break_file,filename,disk,dir,stem,ext
   fname = stem+ext
   year = (byte(strmid(strupcase(fname),0,1))+1929)(0)
   doy = strmid(fname,1,3)

   ; Loop through and analyze all selected scans
   for i = 0, nacal-1 do begin
      ; Check if this scan was selected, and not already analyzed
      if (scan[i].sel eq '*' and scan[i].anal ne '*') then begin
         outfile = !defaults.workdir+string(year,doy,scan[i].srec,format='(I4,A3,"_",I5.5)')+'.sav'

         ; If the keyword UPDATE was set, then if this file already exists, do not
         ; reanalyze it.
         if (keyword_set(update)) then begin
            ; Analyze adds an 'i' to the outfile name, so search for such names
            break_file,outfile,disk,dir,stem,ext
            outfilei = disk+dir+stem+'.sav'
            file = findfile(outfilei,count=nfile)
         endif else nfile=0
         if (nfile eq 0) then $
            analyze,filename,scan[i].srec,120.,[scan[i].srec,scan[i].erec],out=outfile

         scan[i].anal = '*'
      endif
   endfor

   ; Run AMPHIT on each scan/baseline

   files=''
   j = 0
   analidx = where(scan.anal eq '*',nanal)  ; List of indexes where ANAL flag is set
   for i = 0, nacal-1 do begin
      if (scan[i].anal eq '*') then files = [files, $
         !defaults.workdir+string(year,doy,scan[i].srec,format='(I4,A3,"_",I5.5)')+'.sav']
   endfor
   nfiles = n_elements(files)-1

   ; Bail out if no PCAL scans found.
   if (nfiles eq 0) then begin
      if (n_elements(nofile) eq 0) then begin
         ans = widget_message('DAILY: No valid PCAL scans in file '+filename,/ERROR)
      endif else begin
         print,'DAILY: No valid PCAL scans in file '+filename
      endelse
      return,0
   endif else begin
      ; Remove null string in first element of FILES string array
      files = files[1:nfiles]
   endelse

   ; Loop over files
   for j = 0, nfiles-1 do begin
      restore,files[j]

      ; Check if the two 27m antennas have enough good data (no sense continuing
      ; if not)
      nxpect = n_elements(avg[0,*,*])
      good = where(finite(avg[0,*,*]),ngood)
      if (ngood ge nxpect/2) then begin
         ; Determine which baselines are present
         nant = scan[analidx[j]].nant
         CASE nant OF
            5: bl_list = [12,14,15,16,24,25,26]
            6: bl_list = [12,14,15,16,17,24,25,26,27]
            7: bl_list = [12,14,15,16,17,18,24,25,26,27,28]
            else: BEGIN
               print,'DAILY: Unexpected number of antennas (',nant,') in file '+files[j]
               return,0
               END
         ENDCASE

         nbl = n_elements(bl_list)

         ; Loop over baselines
         for i = 0, nbl-1 do begin
            ; Determine average amplitudes and phases for each baseline
            amphit,f,avg,bl_list[i],tavg,header,cfg,out,ndfac=ndfac
            bad = where(finite(reform(out)) eq 0,nbad)
;            if (nbad ne 0) then out[bad] = -99
            scan[analidx[j]].amf[*,*,i] = out   ; Save Amp, dAmp, Pha, dPha
         endfor

         ; Correct the phases for ORDER and BEPHZ so that phases will close in WPHZFIT
         phases = reform(scan[analidx[j]].amf[3,*,0:nbl-1])
         for i = 0, nbl-1 do begin
            ; Determine ORDER and BEPHZ for each baseline, and correct for them
            cal_info,bl_list[i],header,cfg,order,bechan,bephz
            ;phases[*,i] = lobe(order*(phases[*,i]-bephz))
         endfor
         scan[analidx[j]].amf[3,*,0:nbl-1] = phases

         ; Save the center time
         tref = total(tavg)/n_elements(tavg)
         scan[analidx[j]].time = strmid(msec2str(tref),1,8)
      endif else begin
         print,'File ',files[j],' skipped because less than half (',ngood*100/nxpect,'%) of the data are non-NaN'
         analidx[j] = -1   ; Flag this index as not analyzed
         if n_elements(iref) ne 0 then begin
            if iref eq j then begin
               print,'Reference scan has been eliminated.  Please try again with a different "IREF" value.'
               return, 0
            endif         
            if iref gt j then iref = iref-1
         endif
      endelse
   endfor

   good = where(analidx ge 0, nanal)
   if (nanal ne 0) then analidx = analidx[good]
   if (nanal eq 0) then begin
      ; No good scans remaining...
      ; If a filename was specified, we are not in interactive mode, so
      ; exit quietly.
      if (n_elements(nofile) ne 0) then begin
         print,'DAILY: No PCAL scans in file with both 27 m antennas operating.'
         return,0
      endif

      ans = dialog_message(['DAILY: No PCAL scans in file with both 27 m antennas operating.','Please select another file, or Cancel.'],/error)
      free_lun,lun
      goto,ask
   endif

   ; Create header lines for table, depending on if 5 elements or 6.
   CASE nant OF
       5: begin
             lines = ["HRec    Time     p12  poff   p14  poff   p15  poff   p16  poff   p24  poff   p25  poff   p26  poff",$
                      "----- -------- ----------- ----------- ----------- ----------- ----------- ----------- -----------"]
          end
       6: begin
             lines = ["HRec    Time     p12  poff   p14  poff   p15  poff   p16  poff   p17  poff   p24  poff   p25  poff   p26  poff   p27  poff",$
                      "----- -------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- -----------"]
          end
       7: begin
             lines = ["HRec    Time     p12  poff   p14  poff   p15  poff   p16  poff   p17  poff   p18  poff   p24  poff   p25  poff   p26  poff   p27  poff   p28  poff",$
                      "----- -------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- -----------"]
          end
       else:
   ENDCASE

   ; Eliminate all scan entries that were not analyzed
   scan = scan[analidx]

   if (n_elements(refscan) eq 0 and n_elements(iref) eq 0) then begin
      ; Compare scans to find the "best" one to use for reference.  It should have the
      ; best S/N and also not be missing too many frequencies
      iref = best_scan(scan)
   endif

   ; Save a copy of the reference scan info (unless provided in REF keyword)
   if (n_elements(refscan) eq 0) then refscan = reform(scan[iref])

   ; Subtract reference phases from phases in all scans, and send result to
   ; WPHZFIT to find global fit.  Then output the result as text (for now).
   refphz = refscan.amf[3,*,0:nbl-1]
   for i = 0, nanal-1 do begin
      scanout = scan[i].amf[*,*,0:nbl-1]

      scanout[3,*,*] = lobe(scanout[3,*,*] - refphz)

      ; Create WINDAT structure to hold the results of PHZFIT
      f = reform(scanout[0,*,0])
      nf = n_elements(f)
      windat = replicate({  x:  !x, $
                            y:  !y, $
                           ap:  fltarr(nf), $
                          phz:  fltarr(nf),$
                         poff:  0.0, $
                         pslp:  0.0, $
                        sigma:  fltarr(2), $
                           dx:    0, $
                           ix:   -2  },nbl)

      ; Run PHZFIT on each baseline
      if (keyword_set(interactive)) then begin
         wphzfit,info=scan[i].srcname,scanout,sn
         openr,savlun,/get_lun,!defaults.workdir+'save.txt'
         line = ''
         readf,savlun,line
         print,'+'+line+'+'
         if line eq 'USE' then begin
            skip = 'USE'
            dat = fltarr(5)
            for k = 0, nbl-1 do begin
               readf,savlun,dat
               windat[k].pslp = dat[2]
               windat[k].poff = dat[1]
               windat[k].sigma = dat[3:4]
            endfor
         endif else begin
            print,scan[i].srcname, scan[i].date, scan[i].time, scan[i].srec, scan[i].erec, format='(A8," scan at ",A10,A8,1X,I5,":",I5," was skipped by user")'
            skip = 'SKIP'
            ; Turn off scan[i].anal (i.e. remove '*')
            scan[i].anal = ' '
         endelse
         free_lun,savlun
      endif else begin
         for k = 0, nbl-1 do begin
            window,k,/pixmap,xsiz=150,ysiz=150
            transfr ={dx:0, ix:-2, sn:sn}
            ibl = k
            ap  = windat[k].ap
            phz = windat[k].phz
            phzfit,transfr,ap,scanout,ibl,phz,res,sigma
            windat[k].x = !x
            windat[k].y = !y
            windat[k].ap   = ap
            windat[k].phz  = phz
            windat[k].poff = res[0]
            windat[k].pslp = res[1]
            windat[k].sigma = sigma
         endfor
      endelse

      if skip eq 'USE' then begin
         scan[i].poff[0:nbl-1] = pclose(windat.poff,windat.sigma[0],cor=cpoff,0)
         scan[i].pslp[0:nbl-1] = pclose(windat.pslp,windat.sigma[1],cor=cpslp,0,/nolobe)

         ; Choose the POFF (lobe) that is closest to zero
         scan[i].poff = lobe(scan[i].poff,/mid0)

         outstr = ''
         for k = 0, nbl-1 do begin
;            wset, k
;            !x = windat[k].x
;            !y = windat[k].y
;            oplot,[0,f,12],scan[i].poff[k]+scan[i].pslp[k]*[0,f,12],linesty=1
            ; Convert sigma in fit parameters to single sigma at 5 GHz
            sigma = sqrt(windat[k].sigma[0]^2 + (windat[k].sigma[1]*5.0)^2)
            ; Write out 5 GHz phase, derived from phase offset and slope
;            outstr = outstr+string((scan[i].poff[k]+5.0*scan[i].pslp[k])/5.,$
;                       nint(sigma)/5.,format='(F7.1,I4)')
            outstr = outstr+string(scan[i].pslp[k],$
                       nint(scan[i].poff[k]),format='(F7.1,I5)')
         endfor

         ; Write out results as a line
         ; (O_DATFILE is a variable restored from the original ANALYSIS save file.
         print,o_datfile,scan[i].srec,scan[i].time,outstr,format='(A12,2X,I5,1X,A8,A)'
         lines = [lines,string(scan[i].srec,scan[i].time,outstr,format='(I5,1X,A8,A)')]
      endif
   endfor
   
   ; Remove any scans that were 'skipped'
   scanidx = []
   for i = 0, nanal-1 do if scan[i].anal eq '*' then scanidx = [scanidx,i]
   nanal = n_elements(scanidx)
   if nanal eq 0 then begin
      print,'DAILY: Error, no scans analyzed!  Nothing to do.'
      return,0
   endif

   scan = scan[scanidx]

   ; Collapse all scans into a single set of phases for each scan by
   ; adjusting each for its phase slope
   sumx = fltarr(nf,nbl)
   sumy = fltarr(nf,nbl)
   for i = 0, nanal-1 do begin
      ; Set NaNs to zero
      bad = where(finite(scan[i].amf) eq 0,nbad)
      if (nbad ne 0) then scan[i].amf[bad] = 0
      ; Calculate phase fits for each frequency/baseline
      p0 = f#scan[i].pslp + replicate(1,nf)#scan[i].poff
      ; Adjust measured phases by subtracting phase fits
      p = reform(scan[i].amf[3,*,*]) - p0

      ; Convert to X,Y
      sumx = sumx + scan[i].amf[1,*,*]*cos(p*!dtor)
      sumy = sumy + scan[i].amf[1,*,*]*sin(p*!dtor)

      ; Set back to NaN
      if (nbad ne 0) then scan[i].amf[bad] = !values.f_nan
   endfor
   pha = atan(sumy,sumx)/!dtor

   refphz = fltarr(nf,nbl)
   for k = 0, nbl-1 do begin
      ; Extract reference phases from REFCAL
      refphz[*,k] = reform((*refcal.pphase)[2,k,*])
      ; Determine ORDER and BEPHZ for each baseline
      cal_info,bl_list[k],header,cfg,order,bechan,bephz
      ; Correct the phases for ORDER and BEPHZ so that phases will close in WPHZFIT
      ;refphz[*,k] = lobe(order*(refphz[*,k]-bephz))
   endfor

   amf = refscan.amf[*,*,0:nbl-1]
   amf[3,*,*] = lobe(pha - refphz)

   if (keyword_set(interactive)) then begin
      wphzfit,info='final reference calibration',amf,sn
   endif else begin
      for k = 0, nbl-1 do begin
         window,k,/pixmap,xsiz=150,ysiz=150
         transfr ={dx:0, ix:-2, sn:sn}
         ibl = k
         ap  = windat[k].ap
         phz = windat[k].phz
         phzfit,transfr,ap,amf,ibl,phz,res,sigma
         windat[k].x = !x
         windat[k].y = !y
         windat[k].ap   = ap
         windat[k].phz  = phz
         windat[k].poff = res[0]
         windat[k].pslp = res[1]
         windat[k].sigma = sigma
      endfor
      openw,savlun,/get_lun,!defaults.workdir+'save.txt'
      printf,savlun,'USE'
      for i = 0, nbl-1 do begin
         printf,savlun,i,windat[i].poff,windat[i].pslp,windat[i].sigma,format='(I2,4F8.1)'
      endfor
      free_lun,savlun
   endelse

   if (keyword_set(noprinter)) then daily_plot,scan,refscan,o_datfile,lines,amf $
                               else daily_plot,scan,refscan,o_datfile,lines,amf,/printer

   ; The LINES string array must be set back to the values appropriate to the MUX
   ; settings in effect, i.e. inverted for ORDER

   ; I have to think about this!



   ; Write the results to the file, unless overridden
   if (not keyword_set(nowrite)) then result = write_daily(filename,lines,scan)

   print,'DAILY: Done with file '+filename
   return,refscan
end