pro argus_mapfsw,outfile,sc1,sc2,fdnum,ifnum,calsc,nint,nfd

;;
;;Calibrates FSW OTF mapping 
;;
;; INPUT 
;;  outfile = output file name
;;  sc1 = Initial scan leg
;;  sc2 = Final scan leg
;;  fdnum = initial beam number
;;  ifnum = spw
;;  calsc = list of VANE scans
;;  nint = number of integrations in scan leg
;;  nfd = total number of beams available
;;
;;
;; Warning set baseline parameters before running
;; e.g., setregion on input example data
;; and modify nfit parameter below as needed

nn = nint
;;nfit, 5

;;
;;Set baseline parameters for project

;;!g.nregion = 4
;;!g.regions = [11900, 13604, 13843, 14061, 14331, 14414, 14685, 17874]

;;
;;

;;ifnum = 0

;;
;; load data
;; argus beams are saved in pairs of two, so to avoid having to load all 8 files, 
;; we determine the corresponding VEGAS Bank for the given fdnum.

;;nfd=16

if (nfd eq 16) then begin
    if ((fdnum eq 8) or (fdnum eq 10)) then begin 
        beam = 'A'
    endif else if ((fdnum eq 9) or (fdnum eq 11)) then begin
        beam = 'B'
    endif else if ((fdnum eq 0) or (fdnum eq 2)) then begin
        beam = 'C'
    endif else if ((fdnum eq 1) or (fdnum eq 3)) then begin
        beam = 'D'
    endif else if ((fdnum eq 12) or (fdnum eq 14)) then begin
        beam = 'E'
    endif else if ((fdnum eq 13) or (fdnum eq 15)) then begin
        beam = 'F'
    endif else if ((fdnum eq 4) or (fdnum eq 6)) then begin
        beam = 'G'
    endif else if ((fdnum eq 5) or (fdnum eq 7)) then begin
        beam = 'H'
    endif
endif else if (nfd eq 14) then begin
    if ((fdnum eq 6) or (fdnum eq 8)) then begin 
        beam = 'A'
    endif else if ((fdnum eq 7) or (fdnum eq 9)) then begin
        beam = 'B'
    endif else if ((fdnum eq 0) or (fdnum eq 2)) then begin
        beam = 'C'
    endif else if ((fdnum eq 1) or (fdnum eq 3)) then begin
        beam = 'D'
    endif else if ((fdnum eq 10) or (fdnum eq 12)) then begin
        beam = 'E'
    endif else if ((fdnum eq 11) or (fdnum eq 13)) then begin
        beam = 'F'
    endif else if ((fdnum eq 4) or (fdnum eq 5)) then begin
        beam = 'G'
    endif

endif

;baseName = string(session, format='(%"AGBT22B_191_%02d")')
;dataFile = string(baseName, baseName, baseName, beam, format='(%"/home/sdfits/%s/%s.raw.vegas/%s.raw.vegas.%s.fits")')
;; dataFile = string(baseName, baseName, beam, format='(%"./%s.raw.vegas/%s.raw.vegas.%s.fits")')

;print, 'dataFile = ', dataFile

;filein, dataFile

;;
;;
;; Writes output data file for each beam 

OUTfile = outfile ;string(session, sc1, ifnum, fdnum, format='(%"KEEP_s%02d_%i_if%02d_fd%02d.fits")')
fileout, OUTfile

;;
;;avoid showing plots which helps with speed

!g.interactive=0

;;
;;-----------------------------------------------------------------
;;
;; Computing Tcal and Vcal

tbg = 2.725

for c=0, n_elements(calsc) - 1 do begin

    gettp, calsc[c], fdnum=fdnum, ifnum=ifnum, quiet=1

    ;; need to get initial plot or will return error message when freezing without a plot
    
    !g.frozen = 1

    ;;

    twarm = !g.s[0].twarm + 273.15                ;; conversion from C to K
    time = !g.s[0].mjd
    freq = !g.s[0].reference_frequency / 1.e9
    elev = !g.s[0].elevation
    getatmos, el=elev, mjd=time, freq=freq, outvals

    if (n_elements(tau) eq 0) then tau=outvals(0)
    tatm = outvals(2)
    am = 1. / sin(elev * !pi / 180.)
    
    tcal = (tatm -tbg) + (twarm - tatm) * exp(tau * am)
    vcal = median(getdata(0))
    
    if c eq 0 then begin
        tcals = [tcal]
        vcals = [vcal]
    endif else begin
        tcals = [tcals, tcal]
        vcals = [vcals, vcal]
    endelse
endfor

tcal = mean(tcals)
vcal = mean(vcals)

;;
;;---------------------------------------------------------------
;;
;;Get Sample scan for example beam for frequency info
;;  Note: REF (sig_state = 0) and SIG (sig_state = 1)
;; 
;;  calculate frequency switch between reference [0] and signal [1]
;

for s=0, 1 do begin
    
    gettp, sc1, plnum=0, fdnum=fdnum, ifnum=ifnum, sig_state=s, quiet=1
    
    if s eq 0 then obsfreq = !g.s[0].reference_frequency else obsfreq = [obsfreq, !g.s[0].reference_frequency]
endfor

freq = obsfreq[1] / 1.e9
fs = (obsfreq[0] - obsfreq[1]) / !g.s[0].frequency_interval

;;
;;---------------------------------------------------------------------
;;
;; loop over scan legs
;

for ii=sc1,sc2 do begin

    chunk = getchunk(scan=ii, fdnum=fdnum, plnum=0, ifnum=ifnum, count=nchunk)

    ;;
    ;;Averaging all data in scan for average OFF 

    for j =0,nchunk-1 do accum, dc=chunk[j]
    ave
    vec=median(getdata(0))

    ;;
    ;;Compute Tsys* = Tcal / (VANE / OFF -1.0)

    mytsys = tcal / (vcal / vec - 1.)    

    ;;Calibrate every sample in scan via Tsys*(SIG-REF)/REF

    for j =0,nn-1 do begin 
        chunk=getchunk(scan=ii,int=j,fdnum=kk,plnum=0,ifnum=ifnum,count=nchunk)

        ;;nchunk should be 2 for the two frequencies phases of FSW
        ;;note chunk cannot get data selected with sig vs ref but sig=T is first
        ;;then sig=F in data listing (need to match convention above for fs to 
        ;;get frequency scale correct)
        
        ;sig_state=0
        set_data_container, chunk[1], buffer=1

        ;sig_state=1
        set_data_container, chunk[0], buffer=0

        subtract, 0, 1
        divide, 0, 1
        scale, mytsys

        ;;Fold data with copy
        copy, 0, 1
        
        ;;shift it
        gshift, fs

        ;subtract it and divide by 2 (i.e., average the two)
        subtract, 1, 0
        scale,0.5

        ;; update header entries
        !g.s[0].units='Ta*'    
        !g.s[0].tsys=mytsys

        ;;Fit and remove baseline (set before starting with example scan)
        ;bshape
        ;baseline

        ;;Set smoothing for project noting VEGAS/sdfits returns 
        ;;specific channels as NaN which can propagate downstream issues 
        ;;with gridder without any smoothing.
        ;;Also could use "replace" to interpolate across known NaN channels

        ;gsmooth, 3, /decimate
;    
        ;; save spectrum to file
        keep

        ;; write some update in the terminal, every 10th spectrum

        if ((j mod 10) eq 0) then begin
            print,'saving ', ', ', fdnum, ':', ii,' /', sc2, '    ',j, ' /', nn - 1
        endif
    endfor
endfor

;;dummy file to reset fileout parameter
fileout,'JUNK.fits'
        ;
;;defaults
;!g.interactive=1
!g.frozen=0
;!g.has_display=1

return
end

;;;Example gridder command
;;;gbtgridder -c 4170:5680 -a 3 -o s01 --noline --nocont KEEP_s01_*.fits
;;;
;;;Note gbtgridder defaults to Gaussian-kernel that smooths data and
;;;output resolution will be ~2^0.5* Beam_fwhm Grid to nearest pixel or
;;;use a different kernel to preserve spatial resolution.