;
; Script output is a CASA script which can be run to clean and collapse the data. 
; Three external file required in cwd (gaincallist-1mmnoflux/planets.txt/bandpasscalibrators.txt)
; A output directory must be created (format is the input filename without any colons)
;
pro calibratesmaswarm
common global
common data_set
common plo
print,'Starting MIR script...'
pl[*].plot_device='ps'
set_plot,'ps'
argsarray=command_line_args(count=nparams)
path=argsarray[0]
scifile=argsarray[1]
rechunk=argsarray[2]
fullscifile=path+scifile
outdir=repstr(scifile,':','')
spawn,'pwd',pwd
filedir=pwd+"/"+outdir

; These values depend on if been rechunked or not
if rechunk gt '1' then begin
   print, ''
   preavg=round(32/rechunk)
   if preavg lt 1 then preavg=1
   ntrim=round((16384/rechunk)*0.065)
   threshold=6
   print, 'This file has been rechunked by ',rechunk,'. Setting ntrim=',strtrim(ntrim,2),' (6.5% of band) and preavg=',strtrim(preavg,2)
endif else begin
   preavg=32
   ntrim=1065
   threshold=10
   print, 'This file has not been rechunked. Setting ntrim=1065 (6.5% of band) and preavg=32'
endelse


; 1) Open a few scans to determine parameters
print,''
print,'#####################################################################'
print,'###### DETERMINE OBSERVING PARAMETERS FOR ', scifile,' ######'
print,'#####################################################################'
print,''
readdata,dir=fullscifile,int=[1,25],/default
select,/pos,/res

; Check if polarization data
result=uti_distinct(bl[pbf].ipol,npol,/many)
if npol eq 4 then begin
   pol=1
   print,'Polarimetry data...needs a different script'
endif
if npol lt 4 then begin
   pol=0
   print,'Not polarimetry data...can be processed'
endif
if (npol lt 1) then begin
   pol=uk
endif
   
; get info on data
   sources=c.source
   recs=c.rec[uti_distinct(bl.irec,nrec,/many)]
   bands=c.band
   sidebands=c.sb
   baselines=c.blcd
   
; get the numbers for each
   numsources=n_elements(sources)
   numrecs=n_elements(recs)
   numbands=n_elements(bands)
   numsidebands=n_elements(sidebands)
   numbaselines=n_elements(baselines)
   
; Check receivers are unique. Dont want to do reduction twice (e.g. 170611_06:53:43).
   recssort=recs[sort(recs)]
   recsuniq=recssort[uniq(recssort)]
   numrecsuniq=n_elements(recsuniq)
   
; Determine reference antenna
   print, ''
   print, 'Determine reference antenna...'
   print, ''
   integs=''
   refant=0
   refant2=0
   refant3=0
   refantlist=[]
   spawn,'grep "0.00000" '+fullscifile+'/antennas | head -1 | cut -c1-2',pad1
   if pad1 then begin
      reads,string(pad1),refantp1
      s1=strsplit(refantp1,'.',/extract)
      refantp1=strtrim(s1[0],2)
      refantlist=[refantlist,refantp1]
      print,'Antenna',refantp1,' on pad 1'
   endif
   
   spawn,'grep "1\.95" '+fullscifile+'/antennas | grep "7\.4" | cut -c1-2',pad2
   if pad2 then begin
      reads,string(pad2),refantp2
      s2=strsplit(refantp2,'.',/extract)
      refantp2=strtrim(s2[0],2)
      refantlist=[refantlist,refantp2]
      print,'Antenna',refantp2,' on pad 2'
   endif
   
   spawn,'grep "2\.7" '+fullscifile+'/antennas | grep "2\.3" | cut -c1-2',pad3
   if pad3 then begin
      reads,string(pad3),refantp3
      s3=strsplit(refantp3,'.',/extract)
      refantp3=strtrim(s3[0],2)
      refantlist=[refantlist,refantp3]
      print,'Antenna',refantp3,' on pad 3'
   endif
   
   spawn,'grep "5\.0" '+fullscifile+'/antennas | grep "2\.5" | cut -c1-2',pad4
   if pad4 then begin
      reads,string(pad4),refantp4
      s4=strsplit(refantp4,'.',/extract)
      refantp4=strtrim(s4[0],2)
      refantlist=[refantlist,refantp4]
      print,'Antenna',refantp4,' on pad 4'
   endif
   
   spawn,'grep "5\.7" '+fullscifile+'/antennas | grep "1\.89" | cut -c1-2',pad5
   if pad5 then begin
      reads,string(pad5),refantp5
      s5=strsplit(refantp5,'.',/extract)
      refantp5=strtrim(s5[0],2)
      refantlist=[refantlist,refantp5]
      print,'Antenna',refantp5,' on pad 5'
   endif
   
   numintegscompare=0
   foreach element,refantlist do begin
      whatrefant = ' "blcd" like "'+element+'" and "wt" gt "0" '
      result=dat_filter(s_f, whatrefant,/reset)
      if result gt 0 then begin
         integs=uti_distinct(in.int,nint,/many)
         numintegs=n_elements(integs)
         if numintegs ge numintegscompare then begin
            if refant eq 0 then begin
               refant=element
               print,' '
               print, 'Setting 1st reference antenna as ', refant 
               print,' '
            endif
            if refant ne element and refant2 eq 0 then begin   
               refant2=element
               print,' '
               print, 'Setting 2nd reference antenna as ', refant2 
               print,' '
            endif else begin        
               refant3=element
                                ;       print, 'Setting 3rd reference antenna as ', refant3
            endelse
            numintegscompare=numintegs
         endif
      endif
      endforeach
      
      print, ' '
      print,"Receivers: ",recsuniq, numrecs, numrecsuniq
      print," "

      recslist=''

;
;############ START BY RX LOOP
;############ STARTING CALIBRATION - OPEN FULL FILE BY RECEIVER
;

      for v=0, numrecsuniq-1 do begin
         rec=recsuniq[v]
         outnamerx=scifile+'_rx'+rec

; only necessary to make list of recs for casa script
;         recslist=recslist+''''+rec+''','


         print,''
         print,'#############################################################################'
         print,'####### STARTING PIPELINE CALIBRATION FOR ', scifile,' RX=', rec, ' #######'
         print,'#############################################################################'
         print,''
         
; Open file
        print,' '
        print,'readdata,dir=',fullscifile,',rx=',rec,',/default'
        print,' '
        readdata,dir=fullscifile,rx=rec,/default
        endelse
         
        print,''
        print, 'Data is loaded'
         
; flag any data labelled as rx=-1
         result=dat_filter(s_f,'"irec" eq "-1"',/no_notify,/reset)
         if result gt 0 then begin
            print,'Flagging data mislabelled as rx=-1'
            flag,/flag
         endif
         select,/pos,/res
         
;  Flag pointing data
         print,''
         print,'***** Flag any pointing data: ', scifile,' *****'
         print,''
         result=dat_filter(s_f,'"integ" lt "10"',/no_notify,/reset)
         if result gt 0 then begin
            flag,/flag
            print,'Pointing data flagged!'
         endif else print,'No pointing data!'
         print,' '
         select,/pos,/res
         
; get info on data
         sources=c.source
         bands=c.band
         sidebands=c.sb
         baselines=c.blcd
         allintegs=uti_distinct(in.int,nint,/many)
         
; get the numbers for each
         numsources=n_elements(sources)
         numbands=n_elements(bands)
         numsidebands=n_elements(sidebands)
         numbaselines=n_elements(baselines)
         numallintegs=n_elements(allintegs)
         
;calc number to divide integs by to get #time scans for one sb,chunk,bl
         numtotalintegs=numsidebands*numbaselines*numbands
         
         print,''
         print,'..............................................'
         print,'................ DATA SUMMARY ................'
         print,'..............................................'
         print,''
         
         select

         
; 1a) Loading gain calibrator information
         print,''
         print,'1) ***** Loading calibrator information *****'
         print,''
         
         print,'1a) Checking for gain calibrators'
         callist = pwd+"/gaincallist-1mmnoflux"
         openr,/get_lun, unit3, callist
         line=''
         gainsource=''
         calibs=''
         
         while not EOF(unit3) do begin
            readf,unit3,line
            calname=STRLOWCASE(strtrim(line))
            
            for i=1, numsources-1 do begin
               if  STRLOWCASE(strcompress(string(sources[i]),/REMOVE_ALL)) eq strcompress(string(calname),/REMOVE_ALL) then begin
                  match = '"source" eq "'+strtrim(sources[i],2)+'" and "wt" gt "0"'
                                ; print, match
                  result=dat_filter(s_l, match,/no_notify,/reset)
                  if result gt 0 then begin
                     sp[psf].igq=1
                     in[pil].sflux='1'
                     gainsource=[gainsource,strtrim(calname,2)]
                     print, string(sources[i]), ' is a gain calibrator'       
                     calibs = [calibs, sources[i]]
                  endif
               endif
            endfor
         endwhile
         
         
; 1b) Loading passband calibrator information
         print,''
         print,'1b) Checking for bandpass calibrators'
         callist = pwd+'/bandpasscalibrators.txt'
         openr,/get_lun, unit1, callist
         line=''
         passsource=[]
         dopasscal=0
         
         while not EOF(unit1) do begin
            readf,unit1,line
            calname=STRLOWCASE(strtrim(line))
            
            for i=1, numsources-1 do begin
               if  STRLOWCASE(strcompress(string(sources[i]),/REMOVE_ALL)) eq strcompress(string(calname),/REMOVE_ALL) then begin
                  
; then check number of scans
                  match = '"source" eq "'+strtrim(sources[i],2)+'" and "wt" gt "0"'
                                ;      print,match
                  result=dat_filter(s_l, match,/reset,/no_notify)
                  if result gt 0 then begin
                     sumintegs=total(sp[psf].integ)
                     timeonbp=sumintegs/numtotalintegs
                     bptimehrs=timeonbp/3600
                     print,'Time on bandpass calibrator is ',strtrim(string(bptimehrs),2),' hours'
; source has to match and have enough time on it (1/2 hour). Poss check against
; strength of source in the future.
                     if timeonbp gt 1800 then begin
                        print, string(sources[i]), ' is a passband calibrator'         
                        sp[psf].ipq=1
                        passsource=[passsource,strtrim(calname,2)]
                        dopasscal=1
                     endif else begin
                        print,'This bandpass calibrator has insufficient flux: ('+strtrim(sources[i],2)+')'
                     endelse
                  endif
               endif
            endfor
         endwhile
         
         
         if dopasscal eq 0 then begin
            print, 'No bandpass calibrator found or insufficient flux. Bandpass calibration will not be performed'
         endif
         
; 1c) Loading flux calibrator information
         dofluxcal=0
         fluxflag=0
         print,''
         print, '1c) Checking for flux calibrators'
         callist = pwd+'/planets.txt'
         openr,/get_lun, unit2, callist
         line=''
         fluxsource=''
         
         while not EOF(unit2) do begin
            readf,unit2,line
            calname=STRLOWCASE(strtrim(line))
            
            for i=1, numsources-1 do begin
               if  STRLOWCASE(strcompress(string(sources[i]),/REMOVE_ALL)) eq strcompress(string(calname),/REMOVE_ALL) then begin
                  print, string(sources[i]), ' is a flux calibrator'
                  fluxflag=1
                  fluxsource=string(sources[i])
                  calibs = [calibs, sources[i]]
                  dofluxcal=1
               endif
            endfor
         endwhile
         
         if dofluxcal eq 0 then begin
            print, 'No flux calibrator found. Flux calibration will not be performed'
         endif
         
         
; 1d) Get array of science targets
         science=''
         for i=1, numsources-1 do begin
            result=where(sources[i] eq calibs, count)
            if count eq 0 then science=[science,sources[i]]
         endfor
         numscience=n_elements(science)
         
         print,''
         print,'Sources summary'
         print,'Science: ', science
         print,'Calibrator: ', calibs
         print,''
         

; make comma separated list of bands for later
         spws=''
         for i=1, numbands-1 do begin
            whatband = '"band" eq "'+strtrim(bands[i],2)+'" and "wt" gt "0"'
                                ;  print,whatband
            result=dat_filter(s_l, whatband,/reset,/no_notify)
            if result gt 0 then begin
               spws=spws+''''+string(bands[i])+''','
            endif
         endfor
         spwnums=repstr(spws,'s','')
         
         select,/pos,/res
         
; 2) Search and fix spikes
         print,''
         print,'***** 2) Search and fix spikes: ', scifile,' RX=',rec,' ******'
         print,''
         for i=0, n_elements(passsource)-1 do begin
            print,"Using ", passsource[i]
            uti_checkspike,source=passsource[i],threshold=threshold,ntrim=ntrim,/fix
            uti_checkspike,source=passsource[i],threshold=threshold,ntrim=ntrim,/fix
         endfor
         
; 3) Regenerate continuum
         print,''
         print, '***** 3) Regenerate continuum: ', scifile,' RX=',rec,' *****'
         print,''
         result=dat_filter(s_f,'"wt" gt "0"',/no_notify,/reset)
         uti_avgband
         
; 4) Flag high Tsys and apply tsys correction
         print,''
         print, '***** 4) Applying Tsys correction: ', scifile,' RX=',rec,' *****'
         print,''
         select,/pos,/re
         result=dat_filter(s_f, ' "tssb" gt "1800" and "wt" gt "0"',/reset,/no_notify)
         if result  gt 0 then begin
            result=dat_filter(s_f, ' "tssb" gt "1800" and "wt" gt "0"',/reset,/no_notify)
            flag,/flag
            print, 'flagging ',result,' scans with Tsys > 1800K'
         endif
         print,' '
         
         select,/pos,/re
         uti_tsys_fix, high=1500,low=100,loose=50,/refit,/auto
         
         print, '***** Generate Tsys figure *****'
         device,filename=filedir+'/FIGS/'+outnamerx+'_tsys.ps',/color,bits=8,decomposed=0
         !P.MULTI=[0,2,2]
         plot_var, x='el', y='tssb',frame_vars='blcd',color_vars='source',frames_per_page=4
         apply_tsys
         
; Regenerate continuum
         print, '***** Regenerate the continuum channel *****'
         result=dat_filter(s_f,'"wt" gt "0"',/no_notify,/reset)
         uti_avgband
         
; 5) Bandpass calibration (if do passcal=1)
         nch=16384              ; for SWARM
         
         if dopasscal eq 1 then begin
            numbpcals=n_elements(passsource)
            for x=0,numbpcals-1 do begin
               select,source=passsource[x],/pos,/res
               sp[psf].ipq=1
            endfor
            
            select,/pos,/res
           ;# mir_save,'test'+rec+'.mir',/new,/nowait
;            mir_restore,'test'+rec+'.mir',/nowait
            print,''
            print,'***** 5) Bandpass calibration (phase): ', scifile,' RX=',rec,' *****'
            print,''
            print,'** Smoothing: 1   ntrim: ',ntrim,'  preavg: ',preavg,' **'
            print,''
            device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',/color,bits=8
            !P.MULTI=[0,2,2]

            pass_cal,cal_type='pha',preavg=preavg, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='no',/defaults
; Check to make sure it worked (first time round loop)
            grep1=''
            grep2=''
            grep3=''
            spawn,'tail -20 '+filedir+'/MIR_CALIBRATION.log | grep -q ''No convergance'' ; echo $?',grep1
            if grep1 eq 0 then begin
                                ; no it didnt work, try again with preavg2
               print,''
               print, 'Fit attempt was unsuccessful. Trying again...'
               print,''
               refant=refant2
               preavg2=preavg*2
               device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',decomposed=0,/color,bits=8
               !P.MULTI=[0,2,2]
               pass_cal,cal_type='pha',preavg=preavg2, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='no',/defaults
               
                                ;check it worked
               spawn,'tail -20 '+filedir+'/MIR_CALIBRATION.log | grep -q ''No convergance'' ; echo $?',grep2
               if grep2 eq 0 then begin
                                ; no it didn't. Try again with preavg3
                  print,''
                  print, 'Fit attempt was unsuccessful. Trying again...'
                  print,''
                  refant=refant2
                  preavg3=preavg*4
                  device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',decomposed=0,/color,bits=8
                  !P.MULTI=[0,2,2]
                  pass_cal,cal_type='pha',preavg=preavg3, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='no',/defaults

                  
                                ;check to see if attempt 3 worked
                  spawn,'tail -20 '+filedir+'/MIR_CALIBRATION.log | grep -q ''No convergance'' ; echo $?',grep3
                  
                  if grep3 eq 0 then begin
                                ; no it didn't. Give up.
                     print,''
                     print, 'Fit attempt was unsuccessful.'
                     print,''
                     device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',decomposed=0,/color,bits=8
                     !P.MULTI=[0,2,2]
                     pass_cal,cal_type='pha',preavg=preavg3, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='no',/defaults
                  endif
                  if grep3 eq 1 then begin
                                ;yes it did. Applying preavg3 solution.
                     device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',decomposed=0,/color,bits=8
                     !P.MULTI=[0,2,2]
                     pass_cal,cal_type='pha',preavg=preavg3, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='yes',/defaults
                  endif
                  
               endif
               
               if grep2 eq 1 then begin  
                                ;yes it did work second time. Apply solution for preavg2
                  print,''
                  print,'Fit attempt was successful. Running again with solutions applied.'
                  print,''
                  device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',decomposed=0,/color,bits=8
                  !P.MULTI=[0,2,2]
                  pass_cal,cal_type='pha',preavg=preavg2, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='yes',/defaults
               endif
               
            endif
            
            if grep1 eq 1 then begin  
                                ; yes it did work. Apply solution preavg
               print,''
               print,'Fit attempt was successful. Running again with solutions applied.'
               print,''
               device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-pha.ps',decomposed=0,/color,bits=8
               !P.MULTI=[0,2,2]
               pass_cal,cal_type='pha',preavg=preavg, smoothing=1,ntrim=ntrim,tel_bsl='telescope',refant=refant, apply_solution='yes',/defaults
            endif
            
            print,''
            print, 'Regenerate the continuum.....'
            print,''
            result=dat_filter(s_f,'"wt" gt "0"',/no_notify,/reset)
            uti_avgband
            
            
; passcal amp
            print,''
            print,'***** 6) Bandpass calibration (amplitude): ', scifile,' RX=',rec,' *****'
            print,''
            command='"nch" eq "'+strcompress(string(nch),/remove)+'" and "wt" gt "0"'
            print,command
            result=dat_filter(s_f,command,/no_notify,/reset)
;            print,'** Smoothing:',floor(nch/10.)
            device,filename=filedir+"/FIGS/"+outnamerx+'_passcal-amp.ps',decomposed=0,SCALE_FACTOR=1,/color,bits=8
            !P.MULTI=[0,2,2]
            pass_cal,cal_type='amp',smoothing=1,refant=refant,ntrim=ntrim, preavg=preavg,apply_solution='yes',/defaults ;,/noplot
            
            print,''
            print, 'Regenerate the continuum.....'
            print,''
            result=dat_filter(s_f,'"wt" gt "0"',/no_notify,/reset)
            uti_avgband
            
; end if dopasscal=1
         endif else begin
            print, '5-6) Skipping band pass calibration.'
         endelse
         
         
; 8) Phase gain calibration
; 8.1) if there is flux calibrator, do flux cal and measure the fluxes of gain calibrators
         
         if fluxflag eq 1 then begin
            
            print,''
            print,'***** 7) Gain Calibration (phase): ', scifile,' RX=',rec,' *****'
            print,''
            
; for all sources except flux cal do gain_cal without connect (realistic fit)
            result=dat_filter(s_f,'"source" ne "'+fluxsource+'" and "wt" gt "0"',/no_notify,/reset)
            device,filename=filedir+"/FIGS/"+outnamerx+'_gaincal-pha.ps',decomposed=0,/color,bits=8
            !P.MULTI=[0,2,2]
            gain_cal,cal_type='pha',x='hours',smoothing=0.1,tel_bsl='telescope',refant=refant,/defaults ;,/noplot
            
                                ; now do flux cal only with connect (force to 0 - not real but good SNR out)
            result=dat_filter(s_f,'"source" eq "'+fluxsource+'" and "wt" gt "0"',/no_notify,/reset)
            sp[psf].igq=1       ;yes it will be selected as a gain calibrator
            in[pif].sflux=1     ;and its flux is 1
            gain_cal,cal_type='pha',/connect,tel_bsl='telescope',refant=refant,/non_point,/defaults,/noplot
            
; flux measurement
            print,''
            print,'***** 8) Flux Measurement (each sideband will be treated independently): ', scifile,' *****'
            print,''
            
                                ; START PER SIDEBAND LOOP
            for y=0, numsidebands-1 do begin
               print,' ' 
               print, 'Starting ',sidebands[y],'sb'
               
               result=dat_filter(s_f,'"sb" eq "'+sidebands[y]+'" and "band" eq "c1" and "wt" gt "0"',/no_notify,/reset)
 ;;;;;;; result=flux_scale(fluxsource,'c1',flux_inp=flux_inp,/noprint,/scale_bsl)
               result=flux_scale(fluxsource,'c1',flux_inp=flux_inp,/scale_bsl,/noprint)
               result=cal_apply(gain='amp')
               result=sma_flux_measure('c1',names=names,flux=flux,/defaults)
               print,' '
                                ; now select flux cal only and use it to get elev limits
               result=dat_filter(s_f,'"sb" eq "'+sidebands[y]+'" and "source" eq "'+fluxsource+'" and "wt" gt "0"',/no_notify,/reset)
               sp[psf].igq=0
               minel=strcompress(string(floor(min(in[pif].el))),/remove)
               maxel=strcompress(string(ceil(max(in[pif].el))),/remove)
                                ; select all sources again with elev limits
               result=dat_filter(s_f,'"wt" gt "0" and "band" eq "c1" and "el" gt "'+minel+'" and "el" lt "'+maxel+'"',/no_notify,/reset)
               result=sma_flux_measure('c1',names=names2,flux=flux2,snr=snr2,/defaults)
               
               select,/pos,/res,sideband=sidebands[y]
               
               print,' '
               print,'Fluxes summary for sideband=',sidebands[y],'sb...'
               
               for z=1, n_elements(gainsource)-1L do begin
                  iflux=where(names eq gainsource[z])
                  gainflux=flux[iflux]
                  list = '"sb" eq "'+sidebands[y]+'" and "source" eq "'+gainsource[z]+'" and "wt" gt "0"'
                  print,list
                  result=dat_filter(s_f,list,/reset,/no_notify)
                  in[pif].sflux=gainflux[0]
                  print,'The flux for Gain calibrator ',gainsource[z],' is derived as ',gainflux[0],' Jy'
                  
                  iflux=where(names2 eq gainsource[z],count)
                  if count gt 0 then begin
                     gainflux=flux2[iflux]
                     snr=snr2[iflux]
                     if snr gt 20 then begin
                        in[pif].sflux=gainflux[0]
                        print,'The flux for gain calibrator ',gainsource[z],' is derived as ',gainflux[0],' Jy between elevation ', minel, ' and ',maxel
                     endif
                  endif
               endfor
               
               
; 10) Amplitude gain calibration
               
               result=dat_filter(s_f,'"sb" eq "'+sidebands[y]+'" and "wt" gt "0"',/no_notify,/reset)
               print,''
               print,'***** 9) Gain Calibration - amplitude: ', scifile,' ', sidebands[y],'sb *****'
               print,''
               device,filename=filedir+"/FIGS/"+outnamerx+'_'+sidebands[y]+'_gaincal-amp.ps',decomposed=0,/color,bits=8
               !P.MULTI=[0,2,2]
               gain_cal,cal_type='amp',x='hours',smoothing=0.3,/preavg,tel_bsl='telescope',refant=refant,/defaults ;,/noplot
               
            endfor              ;end per sideband loo[
            
         endif else begin
                                ; THIS IS THE ELSE NOT DOING FLUX MEASURE LOOP
            result=dat_filter(s_f,'"wt" gt "0"',/no_notify,/reset)
            print,''
            print,'***** 7) Gain calibration (phase): ', scifile,' *****'
            print,''
            device,filename=filedir+"/FIGS/"+outnamerx+'_gaincal-pha.ps',/color,decomposed=0,bits=8
            !P.MULTI=[0,2,2]
            gain_cal,cal_type='pha',x='hours',smoothing=0.1,tel_bsl='telescope',refant=refant,/defaults ;,/noplot   
            
            print,''
            print,'***** 8-9) Skipping flux measure & gain calibration (amplitude): ', scifile,' RX=',rec,' *****'
            print,''
            
                                ; Dont have proper value for this. Can
                                ; use avg value from bandpass for each
                                ; quasar and use that. Or just decide
                                ; that the Tsys calibration will suffice
;   print,'***** 9) Gain Calibration - amplitude (average flux only - no flux calibration performed): ', scifile,' ', sidebands[y],' *****'
;   print,''
;   device,filename=filedir+"/FIGS/"+outname+'_'+sidebands[y]+'_gaincal-amp.ps',decomposed=0,/color,bits=8
;   !P.MULTI=[0,2,2]
;   gain_cal,cal_type='amp',x='hours',smoothing=0.3,/preavg,tel_bsl='telescope',refant=refant,/defaults ;,/noplot
            
         endelse
         
print, ' '
print,'***** 10) Regenerating the continuum: ', scifile,' RX=',rec,' *****'
print,''
print, 'Regenerate the continuum.....'
print,''
result=dat_filter(s_f,'"wt" gt "0"',/no_notify,/reset)
select,/pos,/re         
uti_avgband


; 10) Combine sidebands and save result
print,''
print,'***** 11) Saving calibrated data for ', scifile,' RX=',rec,' *****'
print,''

filenamecal=filedir+"/CAL/"+outnamerx+'_CAL.mir'
select,/pos,/res
mir_save, filenamecal,/new,/nowait
select,/pos,/res,source='-'+fluxsource
device,filename=filedir+"/FIGS/"+outnamerx+'_cal.ps',decomposed=0,/color,bits=4
!P.MULTI=[0,1,1]
plot_continuum

select,/pos,/res

callist = pwd+"/gaincallist-1mmnoflux"
openr,/get_lun, unit4, callist

for k=1, numscience-1 do begin      
print,'Starting ',science[k]   
   match =''
   spawn,'grep -q '+science[k]+' '+pwd+'/gaincallist-1mmnoflux ; echo $?',match
   print,'match = ',match
 
   if match eq 1 then begin
      autofits,source=science[k]
   endif
endfor                          ; end per science loop

endfor                          ; end per receiver loop
     
endelse                      ; end 'if not pol data loop
print,"end for pol=0"
pl[*].plot_device='x'
print,''
print,'***** Calibration script complete for ', scifile,' *****'
print,''
   
end                             ; end pro