;+
; NAME:
;     ADJUST_AMF
; PURPOSE:
;     Compares multiple .AMF files taken with different baselines
;     during a standard series of reference calibrations, makes
;     the multiple measurements on the 12 baseline consistent,
;     averages the 12 baseline data, and writes out a final .AMF
;     file containing the multiple baseline information.
; CATEGORY:
;     OVRO APC CALIBRATION SUPPORT
; CALLING SEQUENCE:
;     adjust_amf[,files][,outfile=outfile][,/debug]
; INPUTS:
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     files     the names of the input files containing the AMPHIT
;                 data from the series of calibration observations
;                 (output from WAMPHIT).  If omitted, the user is
;                 asked to choose the series of input files.
;     outfile   the name of the output .AMF file to write the
;                 combined information to.  If omitted, the user
;                 is asked to provide a filename interactively.
;     debug     a switch that, when set, causes plots to be shown
;                 for gauging the quality of the result.
; ROUTINES CALLED:
;     lobe
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
;     Writes the results to an output file specified by the user.
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 29-Jul-2000 by Dale E. Gary
;     16-Sep-2000  DG
;       Fixed up several things such as converting -99.9 to NaN on
;       input, and NaN to -99.9 on output.
;-

pro adjust_amf,files,outfile=outfile,debug=debug

   ; See if FILES was set, and if not, ask user for files
   if (not keyword_set(files)) then begin
      files = dialog_pickfile(/multi,path=!defaults.workdir,filter='*.amf',$
                title='Choose multiple files for input (any order)')
      if (files[0] eq '') then return
   endif

   nfiles = n_elements(files)
   nf = (nbl = (ns = 0))  ; Declare type (i.e. integer) of these variables

   ant = intarr(nfiles)
   ; Go through files one at a time
   for i = 0, nfiles-1 do begin
      openr,lun,/get_lun,files[i]
      readf,lun,nf,nbl,ns
      if (ns ne 1) then begin
         print,'ADJUST_AMF: Cannot handle files with more than one scan:'
         print,'            ',files[i]
         free_lun,lun
         return
      endif
      in = {hed:strarr(3),data:fltarr(6,nf)}
      if (i eq 0) then begin
         calin = fltarr(6,nf,nbl)
         calout = fltarr(6,nf,nbl,nfiles)
      endif
      for j = 0, nbl-1 do begin
         readf,lun,in
         calin[*,*,j] = in.data
      endfor
      free_lun,lun

      bad = where(calin eq -99.9,nbad)
      if (nbad ne 0) then calin[bad] = !values.f_nan

      ; Force phase closure of each scan of data
      calout[*,*,*,i] = cal_pclose(calin)  ; This will remove NaN in the phase column

      ant[i] = strmid(in.hed[0],1,1)  ; Read antenna number from last baseline
   endfor

   ; We have read in all of the data.  Now average the 12 baseline
   ; amplitudes and adjust the 12 baseline phases for best agreement
   ; and average them also.

   ; Find the best s/n
   mins = intarr(nf)
   ; Get array of indexes where minimum s/n as a function of frequency
   for i = 0, nf-1 do begin
      junk = min(calout[2,i,0,*]/calout[1,i,0,*],imin)    ; Min s/n vs. freq.
      mins[i] = imin   ; Save index of min s/n
   endfor
   ; Count all of the times index i = 0, nfiles-1 occurs in the mins array
   nmins = intarr(nfiles)
   for i = 0, nfiles-1 do begin
      junk = where(mins eq i,nmin)  ; Where the min occurs for file I
      nmins[i] = nmin               ; Save the number of min for file I
   endfor
   ; Find which index occurs most often
   junk = max(nmins,ibest)

   ; Set up for plotting comparisons, if /debug was used
   if (keyword_set(debug)) then begin
      oldp = !p
      !p.multi = [0,nfiles,2,0,0]
      xy = get_screen_size()
      window,/free,xsiz=xy[0]/1.5,ysiz=xy[1]/(1.5*nfiles/2)
   endif

   f = reform(calout[0,*,0,0])
   ; Use file with index IBEST as the reference for the 12 baseline
   res = fltarr(2,nfiles)
   for i = 0, nfiles-1 do begin
      p = reform(lobe(calout[3,*,0,i]-calout[3,*,0,ibest],/mid0))
      dp = reform(calout[2,*,0,i]/calout[1,*,0,i])
      good = where(finite(p),ngood)
      if (ngood le 3) then begin
         print,'ADJUST_AMF: Not enough points for linear fit, for file ',files[i]
         return
      endif
      res[*,i] = linfit(f[good],p[good],sdev=dp[good],chi=chi,sigma=sigma)
      if (keyword_set(debug)) then begin
         break_file,files[i],disk,dir,stem,ext
         plot,f[good],p[good],psym=1,xtit='Frequency [GHz]',ytit='Phase [degrees]',$
                          tit='Fit BL12 for file '+stem+ext,yran=[-180,180]
         oplot,f,res[0,i]+res[1,i]*f
      endif
   endfor

   ; Apply these slopes to the 12 baseline phases, add half of the slopes
   ; to each of the other baselines to preserve phase closure
   for i = 0, nfiles-1 do begin
      calout[3,*,0,i] = lobe(calout[3,*,0,i] - res[0,i] - res[1,i]*f)
      for j = 1, nbl-1 do begin
         calout[3,*,j,i] = lobe(calout[3,*,j,i] - (res[0,i] + res[1,i]*f)/2)
      endfor
   endfor
   ; Then average the amplitudes and phases through a vector average
   x = reform(calout[1,*,0,*]*cos(calout[3,*,0,*]*!dtor))
   y = reform(calout[1,*,0,*]*sin(calout[3,*,0,*]*!dtor))
   xavg = total(x,2)/nfiles
   yavg = total(y,2)/nfiles
   calout[1,*,0,0] = sqrt(xavg^2 + yavg^2)                   ; Average amplitudes
   calout[3,*,0,0] = lobe(atan(yavg,xavg)/!dtor)             ; Average phases
   calout[2,*,0,0] = sqrt(total(reform(calout[2,*,0,*])^2,2))/nfiles  ; Average rms error

   ; Update phase closure after the above phase adjustments
   for i = 1, nfiles-1 do calout[3,*,0,i] = calout[3,*,0,0]
   for i = 0, nfiles-1 do calout[*,*,*,i] = cal_pclose(calout[*,*,*,i])
   if (keyword_set(debug)) then begin
      for i = 0, nfiles-1 do begin
         break_file,files[i],disk,dir,stem,ext
         plot,f,lobe(calout[3,*,0,i]-calout[3,*,1,i]-calout[3,*,2,i]),psym=1,$
            xtit='Frequency [GHz]',ytit='Phase [degrees]',$
                          tit='Phase closure for file '+stem+ext
      endfor
      !p = oldp
   endif

   ; Convert NaN back to -99.9 so that files can be read with Free-format read
   bad = where(finite(calout) eq 0,nbad)
   if (nbad gt 0) then calout[bad] = -99.9

   ; Ready to write out the results to a new .AMF file
   if (not keyword_set(outfile)) then begin
      outfile = dialog_pickfile(path=!defaults.workdir,filter='*.amf',$
                title='Choose file for output')
      if (outfile eq '') then outfile = !defaults.workdir+'adjust_amf.amf'
   endif

   if (oktowrite(outfile) ne 'Yes') then return

   openw,lun,/get_lun,outfile
   printf,lun,nf,nfiles*2+1,1
   iorder = sort(ant)

   pslope = 0
   pol = 'I'
   ; Print out averaged 12 baseline info
   printf,lun,12,' baseline [',pslope,']',pol,format='(I2,A11,F5.1,A1,1X,A1)'
   printf,lun,' GHz    Amp    dAmp      Pha   dPha   pslope'
   printf,lun,'----  ------- -------  ------ ------  ------'
   for j = 0, nf-1 do begin
      printf,lun,calout[*,j,0,0],format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))'
   endfor

   for i = 0, nfiles-1 do begin
      ; Print out 1n baseline info
      printf,lun,10+ant[iorder[i]],' baseline [',pslope,']',pol,format='(I2,A11,F5.1,A1,1X,A1)'
      printf,lun,' GHz    Amp    dAmp      Pha   dPha   pslope'
      printf,lun,'----  ------- -------  ------ ------  ------'
      for j = 0, nf-1 do begin
         printf,lun,calout[*,j,1,iorder[i]],format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))'
      endfor
   endfor
   for i = 0, nfiles-1 do begin
      ; Print out 2n baseline info
      printf,lun,20+ant[iorder[i]],' baseline [',pslope,']',pol,format='(I2,A11,F5.1,A1,1X,A1)'
      printf,lun,' GHz    Amp    dAmp      Pha   dPha   pslope'
      printf,lun,'----  ------- -------  ------ ------  ------'
      for j = 0, nf-1 do begin
         printf,lun,calout[*,j,2,iorder[i]],format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))'
      endfor
   endfor
   free_lun,lun

return
end