;+
; NAME:
;       calib_rate
; PURPOSE:
;       reads output of analyze & apply calibration and rearrange
;   visibilities into complex array with index (i,j,k) where
;       i = time
;       j = baseline (j = 0,1,2,3,4,5,6,7,8,9)
;       k = frequency
; CATEGORY:
;   OVSA APC imaging
; CALLING SEQUENCE:
;       cal=calib_rate('inname')
;    cf.    calib, 'inname', 'outname'
; INPUTS:
;   infiles <= output of analyze.pro
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;       None at present
; ROUTINES CALLED:
;   calib_ap.pro
;   calib_apply.pro
;   ovsa_poln
;   and some other SSW files
; OUTPUTS:
;       outfile =>  an idlsave file which will input to IMAGR.PRO
; COMMENTS:
;   visibilies are in double flt, but here converted to single flt
;   may need some more actions like 'combine', 'backgr-subtr'
;   before going to IMAGR.PRO
; SIDE EFFECTS: -
; RESTRICTIONS:
;   R and L visibilities upto 7 baselines.
; MODIFICATION HISTORY:
;     Written 05-JAN-2001 by JL
;     20-Mar-2001  DG
;       Changed to use new file convention--all polarizations in the
;       same file.
;     11-Apr-2001  DG
;       Changed to more general form, works with any set of antennas
;       and/or polarizations
;     23-Jul-2003  DG
;       Eliminate redundant code with OVSA_POLN routine, and call that
;       routine.
;     19-Sep-2010  DG
;       Discovered error in calculation of RA and Dec at start time.  RA
;       and Dec are not used in subsequent imaging, but these should be
;       correct since they appear in image FITS headers.
;-

function calib_rate, inname

restore,inname

print,' - apply calibration

avg = apply_cal_all(tavg,avg,f,refcal)

; * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
; Note on Channel numbers vs. Baseline indexes
; Channel numbers are defined as
; BSL:  12   14   15   16    24    25    26     45   46    56
; CHL:  5-6  7-8  9-10 11-12 13-14 15-16 17-18 19-20 21-22 23-24
;     Sin Cos
; POL:  LR   LI   LI   LI    RI    RI    RI     I    I     I
;
; From this point on, we will NOT use channels but only baselines.
; Channel has been used to hold sine & cosine separately
; and one need to figure out which pair of (two adjacent) channels
; corresponds to which baseline.
; This feature is no longer necessary if use complex variables
; or still fltarr with one more index 0=real & 1=imaginary.
; Since we know enough  how IDL does complex var. we use the former.
;
; * * * Index convention * * * *
; time      0 < i < nt
; baseline  0 < j < nb
; frequency 0 < k < nf
;
;* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

nant = refcal.nant
navg = n_elements(avg[0,0,*,0])
nf   = n_elements(avg[0,*,0,0])
npol = n_elements(avg[0,0,0,*])
nbl  = refcal.nbl
ns = nant-2   ; Number of small antennas
t0 = tavg[0]  ; Start time in ms
tavg=tavg/3.6E6  ; Convert times to decimal hours

; Convert data from R, L and I to just R and L (requires that background
; has already been subtracted)
ovsa_poln,avg,nant,nbl,npol,R_,L_,TP_R,TP_L

; Switch order of indexes so that time comes first
R = complexarr(navg,nbl,nf)
L = complexarr(navg,nbl,nf)
TP = complexarr(navg,nant,nf)
for i = 0, navg-1 do begin
   R[i,*,*] = R_[*,*,i]
   L[i,*,*] = L_[*,*,i]
   TP[i,*,*] = complex(TP_R[*,*,i],TP_L[*,*,i])
endfor
R_ = 0  ; Just to free the memory
L_ = 0  ; Just to free the memory
TP_R = 0
TP_L = 0

;-----------------Weighting business (SCBONG...)

sn = (*refcal.pfactors)/(*refcal.prms)
sn_r = reform(sn[0,nant:nant+nbl-1,*])
sn_l = reform(sn[1,nant:nant+nbl-1,*])

WGT = complexarr(navg,nbl,nf)

for i = 0, navg-1 do begin
   WGT[i,*,*] = complex(sn_r^2/reform(R[i,*,*]^2),sn_l^2/reform(L[i,*,*]^2))
endfor

uv = complexarr(navg,nbl,nf)
; for uv we know that there would be no invalid points.
for k = 0, nf-1 do begin
   ibl = -1
   for i=0,nant-2 do begin
      for j=i+1,nant-1 do begin
         ibl = ibl+1
         u= reform(uvt[j,i,*],navg)*f[k]
         v= reform(uvt[i,j,*],navg)*f[k]
         uv[*,ibl,k] = complex(u,v)
      endfor
   endfor
endfor

; other information that needs to be saved
; date='year/mm/dd'
; src='AR0000'
; gtab=[X_cen_sun_arcsec, Y_cen_sun_arcsec, $
;   HA_deg,DEC_deg, PA_deg, B0_deg]
date=tls.date
tstart=strtrim(tls.timstr,2)
src =strcompress(traj.srcname,/remove_all)
tarray  =[1D,t0,t0^2]     ; units of tavg OK?
RA_deg  = total(geometry.a_ra*tarray)/1000./3600.*15.  ; RA at start time (decimal degrees)
DEC_deg = total(geometry.a_dec*tarray)/1000./3600      ; Dec at start time (decimal degrees)
Xcen_as = (*traj.pha)*3.6       ; milli degree -->  arcsec
Ycen_as = (*traj.pdec)*3.6      ; milli degree -->  arcsec
PA_deg  = geometry.solar.pa[0]/10000.  ; milli degree --> degrees
B0_deg  = geometry.solar.b0[0]/10000.  ; milli degree --> degrees
gtab = float([Xcen_as,Ycen_as, PA_deg,B0_deg, RA_deg,DEC_deg])
print, 'Original Datafile = ',O_DATfile
info={src:src, date:date, start:tstart, gtab:gtab, o_datfile:O_DATFILE}
cal={info:info, t:tavg, f:f, uv:uv, RR:R, LL:L, TP:TP, WGT:WGT}
return,cal
end