;+
; NAME:
;     BASEPHZ
; PURPOSE:
;     Allows interactive adjustment of frequency-dependent phase
;     slope and offset in a single "baseline" measurement (a scan
;     processed with ANALYZE, and creates or adds a line to an
;     output file suitable for input to the BASFIT baseline fitting
;     routine.
; CATEGORY:
;     OVSA APC CALIBRATION
; CALLING SEQUENCE:
;     basephz,filename[,calfile=calfile]
; INPUTS:
;     filename   The name of the .SAV file created by ANALYZE, containing
;                  processed (averaged) data from a PCAL scan, organized
;                  into a time series of spectral information.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     calfile    The name of a suitable phase and amplitude calibration
;                  file produced by AMPHIT.  The file contains amplitude,
;                  phase, and uncertainties in amplitude and phase, for
;                  a list of frequencies, for a list of baselines IN ORDER
;                  [12, 14, 15, 16, 24, 25, 26].  If omitted, or included
;                  but equal to an empty string [''], the user
;                  is asked to input a filename.
; ROUTINES CALLED:
;     cal_pclose, amphit, wphzfit
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 01-Jul-1999 by Dale E. Gary
;     31-Mar-2000  DG
;       Significant cleanup and documentation.  Now writes a standard format
;       file that no longer needs editing!
;     14-Apr-2000  DG
;       Get rid of C:\TEMP references (use !Defaults.workdir)
;     10-May-2000  DG
;       Change handling of AMPHIT to reflect change of format of OUT array
;       returned by AMPHIT.
;     05-Jul-2000  DG
;       Now uses a variable NF=6 rather than the number 6, but the number of
;       frequencies is still hardwired, just in only one place now.  Probably
;       the CALFILE (just an AMPHIT file) should contain the number of frequencies.
;     01-Jan-2001  DG
;       Changed to read the number of frequencies, baselines, and scans from
;       the .CAL file (Just an AMPHIT file).  Also, if the .CAL file contains
;       more frequencies than the data, that is okay--only the appropriate
;       frequencies will be used.
;     30-Aug-2004  DG
;       Fixed some things to reflect changes to AMPHIT file (.AMF) format.
;       Also added information string in call to WPHZFIT.
;-
pro basephz,file,calfile=calfile

   ; If the keyword CALFILE is not set, define CALFILE as an empty string
   if (not keyword_set(calfile)) then calfile = ''

   ; If CALFILE is the empty string, ask the user for a CALFILE
   if (calfile eq '') then $
      calfile = dialog_pickfile(filter='*.cal',path=!defaults.workdir,/read)

   ; If still empty, the user must have hit the CANCEL button, so quit.
   if (calfile eq '') then return

   ; Open and read the contents of the CALFILE
   ON_IOERROR,done
   openr,lun,/get_lun,calfile
   readf,lun,nf,nbl,ns
   in = {hed:strarr(4),data:fltarr(6,nf)}
   maxbl = 9
   calin = fltarr(6,nf,nbl)
   bl = intarr(nbl)
   for i = 0, nbl-1 do begin
      readf,lun,in
      bl[i] = fix(strmid(in.hed[0],0,3))
      calin[*,*,i] = in.data
   endfor
done:
   free_lun,lun
   ON_IOERROR,NULL

   ; Restore the input SAVE file specified in the command line
   restore,file

   ; Ensure that the frequencies in the CAL file match those of the data
   for i = 0, nf-1 do begin
      j = where(abs(f - calin[0,i,0]) lt 0.001,nj)
      if (nj eq 0) then calin[*,i,*] = !values.f_nan  ; Mark all data bad for unused freq
   endfor
   nf = n_elements(f)
   ; Eliminate unused frequencies, marked as NaN above.  This will fail if there
   ; are any missing frequencies
   calin = reform(calin[where(finite(calin))],6,nf,nbl)

   ; Open the standard file for writing, appending if necessary
   openw,outlun,/get_lun,!defaults.workdir+'basfit.inp',/append

   ; Check if the file is new, in which case we want to write out
   ; two header lines
   status = fstat(outlun)
   if (status.cur_ptr eq 0) then begin
      ; File is new, so write out the header lines
      lines = ["Filename      HRec   Time   p12  sig   p14  sig   p15  sig   p16  sig   p24  sig   p25  sig   p26  sig",$
               "------------- -----  ----- ---------- ---------- ---------- ---------- ---------- ---------- ----------"]
      printf,outlun,lines,format='(A)'
      print,lines,format='(A)'
   endif


   ; remove any points with zero time (which can occur at the end
   ; of the scans)
   good = where(tavg ne 0,ngood)
   if (ngood gt 0) then begin
      tavg = tavg(good)
      avg = avg[*,*,good]
   endif


   ; Force phase closure of the CALFILE information
   calout = cal_pclose(calin)

   x = intarr(nf)
   junk = ''
   npts = n_elements(tavg)

   ; If the BASELINE scan is extra long, it can be broken into smaller
   ; sections and each used as a BASELINE measurement.  The code below
   ; breaks up scans of more than 19 points into sections of 10 points
   ; each.  A "point" represents data averaged typically over 2 minutes
   ; (the time average is defined in the call to ANALYZE).
   k1 = 0
   k2 = npts-1
   if (npts ge 109) then begin
      ind = indgen(10)*10
      k1 = ind(where(ind lt npts))
      k2 = ((k1+10)<npts)-1
   endif
   ; To avoid excessively short sections, check the last one and if less
   ; than 6 points, combine with the previous section.  If there is only
   ; one section, go ahead with analysis even if it is short.
   nk = n_elements(k1)
   npt = k2(nk-1) - k1(nk-1) + 1  ; Number of points in the last section
   if (nk gt 1 and npt le 5) then begin
      k2(nk-2) = k2(nk-1)
      nk = nk - 1
   endif

;   window,0,xsiz=450,ysiz=600
   blout = calout

   ; Loop over sections in the scan (typically just one section)
   for j = 0, nk-1 do begin
      outstr = ''
      ; Loop over baseline
      for i = 0, nbl-1 do begin
         ; Process each baseline with AMPHIT (just a straight average,
         ; no provision for time drifts).
         amphit,f,avg[*,*,k1[j]:k2[j]],bl[i],tavg[k1[j]:k2[j]],header,cfg,out
         blout(*,*,i) = out
         tim = total(tavg[k1[j]:k2[j]])/(k2[j]-k1[j]+1.)
         tim = tim/1000.  ; Convert time from [msec] to [s]
         ; Apply phase calibration
         blout(3,*,i) = (reform(out(3,*)) - reform(calout(3,*,i))) mod 360
      endfor

      ; WPHZFIT does all the work--displays the data in a widget and allows
      ; the user to adjust the individual phases by +/- 360 degrees to seek
      ; the best global solution, and returns the frequency-dependent phase
      ; information as slope and offset at 1 GHz
      wphzfit,blout,info=strtrim(geometry.srcname)+header.tls.timstr
      openr,tmplun,/get_lun,!defaults.workdir+'save.txt'  ; Read the info written by WPHZFIT
      dat = fltarr(4)
      outstr = ''
      for i = 0, nbl-1 do begin
         readf,tmplun,dat
         ; Write out 5 GHz phase, derived from phase offset and slope
         outstr = outstr+string(dat[1]+5.0*dat[2],nint(dat[3]),format='(F7.1,I4)')
      endfor
      free_lun,tmplun
      ; Organize the data slightly differently and write it out as a line
      ; (O_DATFILE is a variable restored from the original ANALYSIS save file.
      break_file,file,disk,dir,stem,ext
      hrec = strmid(stem,0,5)
      printf,outlun,o_datfile,hrec,long(tim),outstr,format='(A12,2X,A5,I7,A)'
      print,o_datfile,hrec,long(tim),outstr,format='(A12,2X,A5,I7,A)'
   endfor
   free_lun,outlun
return
end