;+
; NAME:
;     ADJUST_REFCAL
; PURPOSE:
;     Applies phase variations at the common frequency (3.6 GHz) to
;     other frequencies, assuming a linear frequency dependence, and
;     writes out a final .AMF file
; CATEGORY:
;     OVRO APC CALIBRATION SUPPORT
; CALLING SEQUENCE:
;     adjust_refcal[,filename]
; INPUTS:
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     filename  the name of the input file containing the AMPHIT
;                 data from the 5-scan calibration (output from
;                 WAMPHIT).  If omitted, the user is asked to choose
;                 an input file.
; ROUTINES CALLED:
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 14-May-2000 by Dale Gary
;     29-Dec-2000  DG
;       Changed to read the number of frequencies, baselines, and
;       scans directly from the .AMF file.
;-
pro adjust_refcal,filename

   if (n_elements(file) eq 0) then begin
      filename = pickfile(path=!defaults.workdir,filter='*.amf')
      if (filename eq '') then return
   endif

   nf = 10    ; Number of frequencies in each scan
   nb = 7     ; Number of baselines (for 5 ants)
   ns = 5     ; Number of scans in REFCAL measurement

   ; Open the file and read through it just to determine what frequencies
   ; are included
   f = 0.
   openr,lun,/get_lun,filename

   readf,lun,nf,nb,ns
   ; Structure of data in an AMPHIT file
   in = {hed:strarr(3),data:fltarr(6,nf)}

   ; Loop over scans
   for j = 0, ns-1 do begin
      for i = 0, nb-1 do begin
         readf,lun,in
         if (i eq 0) then f = [f,reform(in.data[0,*,0])]
      endfor
   endfor
   free_lun,lun
   ; Remove extraneous 0.
   f = f[1:*]
   ; Sort the frequencies
   f = f(sort(f))
   ; Keep only unique frequencies
   f = f(uniq(f))
   ; Total number of frequencies in a REFCAL measurement
   nftot = n_elements(f)

   ; Reopen the file to read the data
   openr,lun,/get_lun,filename

   readf,lun,nf,nb,ns
   ; Array to hold the accumulated output
   refcalout = fltarr(6,nftot,nb)
   cal36 = fltarr(6,nb,ns)  ; To hold 3.6 GHz data

   ; Structure of data in an AMPHIT file
   in = {hed:strarr(3),data:fltarr(6,nf)}
   calin = fltarr(6,nf,nb)
   hedin = strarr(3,nb)

   ; Phase slope corrections as a function of baseline and scan
   pcor = [[  0,  0, 19, 28,  0, 19, 28],$
           [  5,  4, 26, 42, -1, 21, 37],$
           [ 17,  7, 37, 54,-10, 20, 37],$
           [ 22, 29, 59, 80,  7, 37, 58],$
           [-22,-38, -8, 26,-16, 14, 48]]

   for i = 0, ns-1 do begin
      ; Subtract avg of above and convert from 5 to 1 GHz
      pcor[*,i] = (pcor[*,i]-total(pcor,2)/ns)/5.0
   endfor

   ; Loop over scans
   for j = 0, ns-1 do begin
      for i = 0, nb-1 do begin
         readf,lun,in
         calin(*,*,i) = in.data
         if (j eq 0) then begin
            hedin(*,i) = in.hed
         endif
      endfor

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

      ; Apply phase slope corrections, if any
      if (n_elements(pcor) ne 0) then begin
         padj = reform(calin[0,*,0])#reform(pcor[*,j])
         calin[3,*,*] = calin[3,*,*] - padj
      endif

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

      ; Determine the phase difference of 3.6 GHz phases wrt first scan
      if (j eq 0) then begin
         pdiff = fltarr(nb)
         p36 = reform(calout[3,2,*])
      endif else begin
         pdiff = lobe(reform(calout[3,2,*] - p36),/mid)
      endelse
      print,pdiff

      ; Skip adjustment wrt 3.6 GHz if the phase slope correction has been applied
      if (n_elements(pcor) ne 0) then pdiff = pdiff*0

      ; Make an index into the REFCALOUT array for each unique frequency in CALOUT
      ifrq = intarr(nf)
      for k = 0, nf-1 do begin
         ifrq[k] = where(f eq calout[0,k,0])
         ; Adjust phases at each frequency according to 3.6 GHz phase difference
         calout[3,k,*] = (calout[3,k,*] - pdiff*calout[0,k,*]/3.6 $
                 + 18000.) mod 360.
      endfor
      ; Save 3.6 GHz (index 2) data since we have to combine it properly
      cal36[*,*,j] = reform(calout[*,2,*])

      ; Insert data with adjusted phases into appropriate part of REFCALOUT array,
      ; using IFRQ array to place them in the correct slots
      refcalout[*,ifrq,*] = calout
   endfor

   free_lun,lun

   ; Make 3.6 GHz data consistent
   for i = 0, nb-1 do begin
      amp = total(reform(cal36[1,i,*]))/ns
      damp = sqrt(total(reform(cal36[2,i,*])^2))/ns
      dpha = (damp/amp)/!dtor
      kfrq = where(f eq 3.6)
      refcalout[1,kfrq,i] = amp
      refcalout[2,kfrq,i] = damp
      refcalout[4,kfrq,i] = dpha
   endfor

   ; Need to do something with the result, refcalout.
   ; For now, write out a single .AMF file with the 7 baselines
   ; listed for all frequencies

   nf = n_elements(refcalout[0,*,0])
   filename = pickfile(path=!defaults.workdir,filter='*.amf',/write)
   if (filename eq '') then return
   if (oktowrite(filename) eq 'Yes') then begin
      openw,lun,/get_lun,filename
      printf,lun,nf,nb,1
      bad = where(finite(refcalout) eq 0,nbad)
      if (nbad ne 0) then refcalout[bad] = -99.9
      for i = 0, nb-1 do begin
         printf,lun,hedin(*,i),format='(a)'
         for ifrq = 0, nf-1 do printf,lun,refcalout(*,ifrq,i),format='(f4.1,1x,2f8.3,1x,3(1x,f6.1))'
      endfor
      free_lun,lun
   endif

return
end