; to convert from SDFITS file for project to 
; ascii data for PSI-Plot
; 
;------------------------------------------------------------------
; enhancements for doing scan-by-scan calibration
; output file includes "date-obs"
; 
; June 2008 : comparisons with glish maps: readglish, glishcomp,
;             glishcomp2, glishcomp3, glishcomp4
; Aug 2008: add listprocs, changes to writedata for AVpol=1
;           More work on plotcont
; Sep 08 : remove all that glish comparison stuff
; 
; Version of August 2011
;  keeps LST column in cal data output
;  unwraps 360 to 0 degree transitions
;  option /diff in getdata subtracts tracking center (good for MOON).
;---------------------------------------------------------------------

; First, run sdfits on the desired datast
;  e.g.: sdfits -backends=dcr -scans=1:24 /home/gbtdata/AGBT08A_014_03 rosette.fits
; 
; Usage:
;    idl
;    > .com topsi
;    > topsi,"filename"      [select sdfits file]
;    > listprocs
;    > getdata,1,45   [give it the scan number range]
;
;  or if you want calibrated data, then:
;    > getdata,1,45,/cal
;
;  to use all scans in the sdfits file:
;    > getdata,0,0 ,/cal [,/ave, /diffpos ]
;
;       /cal causes output to be in units of Kelvin
;       /ave averages the 2 polarizations.
;       /diffpos : positions relative to a possibly moving center (e.g., the moon)
;
;  alternately:
;   > getscans,"filename",scan1,scan2   : does "topsi" + "getdata,,/cal
;       reads and calibrates specified scans from filename

;
;  to write out the data:
;     > writedata, [ /raw ]
;     [the ascii dump goes into file "outcaldata", or "outrawdata" ]

;  various plots:
;     plotcover   - plots the scan pattern
;     plotdata    - plots all data vs LST
;     plotcont    - tries to do a contour plot
;     gridimage   - tries to grid the image and plot as color scale

; --------------------------------------------------------
; a "typical" session
;  
;  sdfits -backends=dcr -scans=400:406,454:456,510:524 /home/gbtdata/AGBT08A_014_03
;  mv AGBT08A_014_03.raw.dcr.fits SUN.raw.dcr.fits

;  idl
;  > .com topsi
;  > topsi, 'SUN.raw.dcr.fits'
;  > getdata,0,0,/cal
;
;  > plotdata      ; plot the data in time sequence vs LST
;  > plotcont      ; plot data as a contour map
;  > plotcover     ; plot RA vs DEC

;  > writedata     ; write out to an ascii file
;
pro helpinfo
  print, "First, run sdfits on the desired dataset "
  print, "e.g.: sdfits -backends=dcr -scans=1:24 /home/gbtdata/AGBT08A_014_03 rosette.fits"
  print, " "
  print, "idl"
  print, "> .com topsi"
  print, "> topsi,'rosette.fits'      [select sdfits file]"
  print, "> listprocs  -- summarize the scans"
  print, "> getdata,1,24,/cal [,/ave] "
  print, "      /cal calibrates to Kelvin; /ave averages the 2 pols."
  print, " "
  print, "> getdata,0,0 ,/cal   -- use all scans in the sdfits file"
  print, " "
  print, "> writedata, [ /raw ]   -- to write out the data"
  print, '[the ascii dump goes into file "outcaldata", or "outrawdata" ]'
  print, ' '
  print, 'Various plots: '
  print, '   plotcover - plot scan pattern'
  print, '   plotdata  - plot all data vs LST'
  print, '   plotcont  - tries to make a contour diagram'
  print, '   gridimage - tries to interpolate to a grid and plot as a image'
end

;-------------------------------------------------------------
; 
;-------------------------------------------------------------
; filein lets you pick a raw dcr fits file from /home/sdfits/
;
function filein,ffile
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans

  if not keyword_set(ffile) then begin
    ffile = ''
    ffile = dialog_pickfile(path='/home/sdfits', filter='AGBT*.raw.dcr.fits')
  endif

  print, 'Reading FITS file:', ffile
  ff = ffile
  return, mrdfits(ffile,1,head)
end

function date2jd,year,mon,day

   jd = (1461L*(year+4800L+(mon-14L)/12))/4
   jd = jd + (367*(mon-2-12*((mon-14)/12)))/12
   jd = jd - (3*((year+4900+(mon-14)/12)/100))/4
   jd = jd + day - 32075

   ; if( mjd!= (int *)0)  *mjd = jd - 2400001;

   return,jd
end
;
;-------------------------------------------------------------
; getdata gets the RA,DEC,DATA, in the scan range sc1-sc2
; if sc1=sc2=0 then select all data in the input SDFITS file.

; selected data is returned in common "s"
; if /cal, then scan-by-scan temp calibration is done.
; if /ave, then average the 2 polarizations.
; if /base, then subtract a baseline from each scan.
;  diffpos : subtract TRGTLONG and TRGTLAT i.e., positions relative to a
;     (possibly moving) target position.
;
pro getdata,sc1,sc2, cal=cal,ave=ave,base=base, diffpos=diffpos
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff,xonn,yoff,yonn, rr, objct,freqq, fwhm

  if not keyword_set(cal) then cal=0
  if not keyword_set(ave) then ave=0
  if not keyword_set(base) then base=0
  if not keyword_set(diffpos) then diffpos=0

  print,"getdata", cal, sc1, sc2

  ; beam size
  ; obsfreq = d[0].obsfreq   ; this is in Hz
  ; fwhm = 1.2* 206265.0*(3.0e8)/(60.0*100.0*obsfreq)  ; in arcmin
  ; print, "obsfreq=", obsfreq/1.0e9, "  fwhm=",fwhm, " arcmin"

  AVpol=ave

  if sc2 le 0 then sc2 = 100000000
  cname1 = d[0].ctype2
  cname2 = d[0].ctype3

  sizd = size(d)
  nrows = sizd[1]

  ; create arrays for the data
  c1 = dblarr(nrows)     ; coordinates
  c2 = dblarr(nrows)
  dd = strarr(nrows)     ; UT time string
  tt = dblarr(nrows)     ; LST time stamps
  jds = dblarr(nrows)     ; MJDsecs time stamps
  ss = intarr(nrows)     ; scan numbers
  scans = intarr(nrows)  ; list of unique scan nums
  xoff = dblarr(nrows)   ; X-pol Cal off data
  xonn = dblarr(nrows)   ; X-pol Cal on  data
  yoff = dblarr(nrows)   ; Y-pol Cal off data
  yonn = dblarr(nrows)   ; Y-pol Cal on  data

  jj = 0L                 ; index for X-pol
  kk = 0L                 ; index for Y-pol

  iscan = 0L   ; counts unique scan numbers
  oldscan = -1L

  for i=0L,nrows-1,2  do begin

    if d[i].scan ge sc1 and d[i].scan le sc2 then begin

      if d[i].scan ne oldscan then begin
        scans[iscan]=d[i].scan
        iscan=iscan+1
      endif

      oldscan = d[i].scan

      objct = d[i].object
      freqq = d[i].obsfreq

      ; select X or LCP polarization, i.e., crval4=-5 or -2
      if d[i].crval4 eq -5 or d[i].crval4 eq -2 then begin
	if d[i].cal eq 'F' then begin
	   xoff[jj] = d[i].data
	   xonn[jj] = d[i+1].data
	endif else begin
	   xonn[jj] = d[i].data
	   xoff[jj] = d[i+1].data
	endelse

        tt[jj] = d[i].lst  
        dd[jj] = d[i].date_obs
        reads,d[i].date_obs, yr,mon,day,hhr,mmin,ssec, $
		format='(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f5.2)'
        mjd = date2jd(yr,mon,day) - 2400001L
        deltaday = (hhr + ((mmin + (ssec/60.0d))/60.0d) )/24.0d
        ; print, hhr,mmin,ssec,deltaday
        mjd1 = mjd + deltaday
        jds[jj] = mjd1 * 24.0d * 3600.0d
        ; jds[jj] = mjd1 * 24.0d * 3600.0d + d[i].duration
        
        c1[jj] = d[i].crval2
        c2[jj] = d[i].crval3

        ; option: get coordinates relative to center of moving target:
        if diffpos eq 1  then begin
          c1[jj] = d[i].crval2 - d[i].trgtlong
          c2[jj] = d[i].crval3 - d[i].trgtlat
        endif

	ss[jj] = d[i].scan
	TCal1  = d[i].tcal     ; all these tcals should be the same
        ; get the names of the coordinates
	jj = jj+1
      endif

      ; select Y or RCP polarization, i.e., crval4=-6 or -1
      if d[i].crval4 eq -6 or d[i].crval4 eq -1 then begin
	if d[i].cal eq 'F' then begin
	   yoff[kk] = d[i].data
	   yonn[kk] = d[i+1].data
	endif else begin
	   yonn[kk] = d[i].data
	   yoff[kk] = d[i+1].data
	endelse
	TCal2  = d[i].tcal
	kk = kk+1
      endif
    endif
  
  endfor

  ; trim data arrays
  nnrows = min(jj,kk)
  c1 = c1[0:nnrows-1]
  c2 = c2[0:nnrows-1]
  dd = dd[0:nnrows-1]
  tt = tt[0:nnrows-1]
  ss = ss[0:nnrows-1]
  jds = jds[0:nnrows-1]
  xoff = xoff[0:nnrows-1]
  xonn = xonn[0:nnrows-1]
  yoff = yoff[0:nnrows-1]
  yonn = yonn[0:nnrows-1]

  scans = scans[0:iscan-1]

  ; adjust if near 360 deg
  if max(c1)-min(c1) gt 300 then begin
     wzz = where(c1 gt 300)
     c1[wzz] = c1[wzz]-360
  endif

  obsfreq = freqq   ; this is in Hz
  fwhm = 1.2* 206265.0*(3.0e8)/(60.0*100.0*obsfreq)  ; in arcmin
  print, "obsfreq=", obsfreq/1.0e9, "  fwhm=",fwhm, " arcmin"

  ; if cal == 0 then just fill the "xcal" and "ycal" 
  xcal=xoff
  ycal=yoff
  
  ; if cal ==1, then do the calibration.
  if cal eq 1 then begin
    print,'call docal', n_elements(xoff),n_elements(xonn),n_elements(yoff),n_elements(yonn),base
    docal, xoff, xonn, yoff, yonn, base=base
  endif

end



;---------------------------------------------------------
; plot a contour diagram of the selected calibrated data
; pol selects pol 1 or 2
pro plotcont,mincont=mincont, hcopy=hcopy, pol=pol
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff,xonn,yoff,yonn, rr, objct,freqq, fwhm

  if not keyword_set(mincont) then mincont=0
  if not keyword_set(hcopy) then hcopy=0
  if not keyword_set(pol) then pol=0

  titl = string( objct, freqq*1.0e-9, 'GHz', format='(a,1x,f6.2,1x,a)')

  if hcopy eq 1 then begin
    set_plot,'ps'
    device, filename='contourplot.ps', /landscape
  endif

  if pol eq 0 then aa = 0.5*(xcal+ycal)
  if pol eq 1 then aa = xcal
  if pol eq 2 then aa = ycal

  levs = [ 0.14, 0.2, 0.3, 0.4, 0.6, 0.8, 1.2, 1.7, 2.4, 3.4, 4.8, 6.8, 9.6]
  rrange = max(aa) - min(aa)
  llevs = min(aa) + levs*(rrange/10)
  contour, aa, c1, c2, /irregular, levels=llevs, xtitle=cname1, $
	ytitle=cname2, /isotropic, title=titl
  print, "levels=", llevs, " K"

  if hcopy eq 1 then begin
    device, /close
    set_plot, 'x'
    spawn, 'lp -d net contourplot.ps'
  endif
end


; plotdata plots all the data in xcal, ycal
;
pro plotdata, pol=pol, hcopy=hcopy
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans

  if not keyword_set(pol) then pol=0
  if not keyword_set(hcopy) then hcopy=0
  if hcopy eq 1 then begin
    set_plot,'ps'
    device, filename='dataplot.ps', /landscape
  endif

  if pol eq 0 then begin
    plot, tt, xcal, title=ff, xtitle="LST(seconds)", ytitle="Tant(K)"
    oplot, tt, ycal
  endif else begin
    if pol eq 1 then begin
      plot, tt, xcal, title=ff, xtitle="LST(seconds)", ytitle="Tant(K)"
    endif else begin
      plot, tt, ycal, title=ff, xtitle="LST(seconds)", ytitle="Tant(K)"
    endelse
  endelse

  if hcopy eq 1 then begin
    device, /close
    set_plot, 'x'
    spawn, 'lp -d net dataplot.ps'
  endif
end


;---------------------------------------------------------
pro plotcover
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff,xonn,yoff,yonn, rr, objct,freqq, fwhm

  ; ttl=string("coverage, ", ff)
  ttl = "coverage"
  plot, c1,c2, xtitle=cname1, ytitle=cname2, title=ttl, /ynozero

end

;---------------------------------------------------------
;
; list object, obsmode, etc for a range of scan numbers
;
pro list,sc1,sc2
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
  if not keyword_set(sc1) then sc1 = 1
  if not keyword_set(sc2) then sc2 = sc1+20
  lastscan=0
  sizd = size(d)
  nrows = sizd[1]
  firstscan = d[0].scan

  for i=0L,nrows-1 do begin
    if d[i].scan ne lastscan then begin
      if d[i].scan ge sc1 and d[i].scan le sc2 then begin
        print, d[i].scan, d[i].object, d[i].obsmode, d[i].observer, d[i].procseqn, $
		"of", d[i].procsize, $
		format='(i5,1x,a16,1x,a18,1x,a12,1x,i4,1x,a3,1x,i4)'
      endif
    endif
    lastscan = d[i].scan
  endfor
  print, "Scans:",  firstscan, " to ", lastscan
end

; finds contiguous groups of scans
pro listprocs
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
  lastscan=0
  sizd = size(d)
  nrows = sizd[1]
  firstscan = d[0].scan
  lastscan  = 0
  lastobj   = ''
  lastproc  = ''

  print,"..scans..", "source", "procedure", $
     format='(2x, a10, 2x, a16, 1x, a24)'


  for i=0L,nrows-1 do begin
    if d[i].scan ne lastscan then begin
      objct = strcompress(d[i].object, /rem)
      obsproc = strcompress(d[i].obsmode, /rem)
      if objct ne lastobj or obsproc ne lastproc then begin
         if lastscan gt 0 then begin
            print, firstscan, lastscan, lastobj, lastproc, $
		format='(i6, 1x, i6, 1x, a16, 1x, a24)'
            firstscan = d[i].scan
         endif
      endif
    endif
    lastscan = d[i].scan
    lastobj = objct
    lastproc = obsproc
  endfor
  if lastscan gt 0 then begin
    print, firstscan, lastscan, lastobj, lastproc, $
		format='(i6, 1x, i6, 1x, a16, 1x, a24)'
  endif

end

;
;---------------------------------------------------------
; topsi is the main program -- start with this.
pro topsi,ffile
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common xx2, docplots

  ; initialize calibration params
  TCal1 = 0
  TCal2 = 0
  AVpol = 0
  qgalac = 0
  docplots='off'

  d = filein(ffile)
  sizd = size(d)
  nrows = sizd[1]
  firstscan = d[0].scan
  lastscan  = d[nrows-1].scan
  print, nrows, " points read"

  listprocs

  ; getdata,0,0,/cal

end


pro getscans,filename, scan1,scan2
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common xx2, docplots

  topsi,filename
  getdata,scan1,scan2,/cal
end


;---------------------------------------------------------
; write the data to text file
;   pol = select just one pol - zero out the other
;---------------------------------------------------------
pro writedata,raw=raw,secs=secs, pol=pol
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff,xonn,yoff,yonn, rr, objct,freqq, fwhm

  if not keyword_set(raw) then raw=0
  if not keyword_set(pol) then pol=0
  if not keyword_set(secs) then secs=0
  if not keyword_set(deci) then deci=0

  nnrows = n_elements(xcal)
  if pol eq 1 then ycal[0:nnrows-1]=0.0
  if pol eq 2 then xcal[0:nnrows-1]=0.0

  ; print raw data
  if raw eq 1 then begin

     openw,lun, 'outrawdata', /get_lun
     ; print headers
     printf,lun, "## RAW DATA in Counts, ", ff
     printf,lun,cname1,cname2,'SCAN', 'LST', 'XCALOFF', 'XCALON', $
	'YCALOFF', 'YCALON', $
	format='(1x,a9,2x,a9,2x,a6,1x,a9,2x,a10,2x,a10,2x,a10,2x,a10)'

     nnrows = n_elements(xcal)
     for i=0L, nnrows-1 do begin
	printf,lun, c1[i],c2[i],ss[i], tt[i], xoff[i],xonn[i],yoff[i],$
	yonn[i], format='(1x,f9.4,2x,f9.4,2x,i6,1x,f9.2,2x,f10.2,2x,f10.2,2x,f10.2,2x,f10.2 )'
     endfor
     close,lun
     free_lun, lun

  endif else begin

  ; here write the maybe calibrated data
	; writecal,  lun, nnrows,c1,c2,ss, cname1, cname2
    openw,lun, 'outcaldata', /get_lun
    nnrows = n_elements(xcal)
    ; print header
    printf,lun, "## filein=", ff
    printf, lun, "##  Tcals=", TCal1,TCal2, " K"
    printf, lun, "##  Source=", objct
    printf, lun, "##  Date-time=", dd[0]
    printf, lun, "##  Obs Freq=", 1.0e-6 * freqq,  " MHz"
    if AVpol eq 1 then begin                   ; simple file for AVpol
       printf,lun,cname1,cname2, 'SCAN', 'LST', 'DATA',  $
	format='(1x, a9, 2x, a9, 2x, a6, 2x,a9, 2x, a12 )'

       for i=0L, nnrows-1 do begin
         printf,lun, c1[i],c2[i], ss[i], tt[i], xcal[i], $
           format='(1x, f9.4, 2x, f9.4, 2x, i6, 2x, f9.2, 2x, f12.5 )'
       endfor
    
    endif else begin
      if secs eq 0 then begin
       printf,lun,cname1,cname2, 'SCAN', 'LST', 'XDAT', 'YDAT',  $
	format='(1x, a9, 2x, a9, 2x,a6, 2x, a9, 2x, a12, 2x, a12 )'

        for i=0L, nnrows-1 do begin
 	 printf,lun, c1[i],c2[i],ss[i], tt[i], xcal[i],ycal[i], $
	 format='(1x, f9.4, 2x, f9.4,2x,i6, 2x,f9.2, 2x, f12.5, 2x, f12.5 )'
        endfor
       endif else begin
       ; write it with the JD seconds
         printf,lun,cname1,cname2, 'SCAN', 'JDSECS', 'XDAT', 'YDAT',  $
	  format='(1x, a9, 2x, a9, 2x,a6, 1x,a15, 2x, a12, 2x, a12 )'

         for i=0L, nnrows-1 do begin
 	    printf,lun, c1[i],c2[i],ss[i], jds[i], xcal[i],ycal[i], $
	    format='(1x, f9.4, 2x, f9.4,2x,i6,1x, f15.2, 2x, f12.5, 2x, f12.5 )'
         endfor
       endelse
      endelse

    close,lun
    free_lun, lun

  endelse
  
end


;---------------------------------------------------------
; set = 'on' to plot raw data and delta cals; set = 'cal' to plot calibrated data.
pro docalplots,onoroff
common xx2, docplots

  docplots='off'
  if onoroff eq 'on' then docplots='on'
  if onoroff eq 'cal' then docplots='cal'
end

;---------------------------------------------------------
; docal does the simple gain calibration
; scan by scan
pro docal, xoff,xonn,yoff,yonn, base=base
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff_1,xonn_1,yoff_1,yonn_1, rr, objct,freqq, fwhm
common xx2, docplots

  print, "Calibrating, tcals=", TCal1, TCal2
  if not keyword_set(base) then base=0

  ; loadct,12

  ; get median of cal deflection, scan by scan

  nnscans = n_elements(scans)

  sm1pos = fltarr(nnscans)
  sm2pos = fltarr(nnscans)
  c1sdif = 0.0 & c2sdif=0.0
  cd1smdif = 0.0 & cd2smdif=0.0
  
  for i=0,nnscans-1 do begin
    wsc = where( ss eq scans[i])   ; pick out this scan number
    nwsc = n_elements(wsc)
    c1d = c1[wsc]  & c1mean = mean(c1d) & sm1pos[i]=c1mean
    c2d = c2[wsc]  & c2mean = mean(c2d) & sm2pos[i]=c2mean
    c1delt = 60.0*(c1d[nwsc-1] - c1d[0])/nwsc  ; integration size in arcmin
    c2delt = 60.0*(c2d[nwsc-1] - c2d[0])/nwsc
    cd1smdif = cd1smdif+abs(c1delt)  & cd2smdif=cd2smdif+abs(c2delt)

    if i gt 0 then begin           ; calc mean coord diff between scan centers, arcmin
      c1sdif = 60.0*(sm1pos[i] - sm1pos[i-1])
      c2sdif = 60.0*(sm2pos[i] - sm2pos[i-1])
    endif

    xdifz = (xonn[wsc]-xoff[wsc])
    xdif = median( (xonn[wsc]-xoff[wsc]))
    ; xdif = mean( (xonn[wsc]-xoff[wsc]))
    xcal[wsc] = (TCal1/xdif) * 0.5*( (xonn[wsc]+xoff[wsc]) - xdif)

    ydifz = (yonn[wsc]-yoff[wsc])
    ydif = median( (yonn[wsc]-yoff[wsc]))
    ; ydif = mean( (yonn[wsc]-yoff[wsc]))
    ycal[wsc] = (TCal2/ydif) * 0.5*( (yonn[wsc]+yoff[wsc]) - ydif)

    print, "Cal ", scans[i], xdif, ydif, nwsc, c1delt,c2delt, c1sdif, c2sdif, $
	format='(a, 2x, i4, 2x, f8.2, 1x, f8.2, 2x, i5, 2x, f8.4, 2x, f8.4, 2x,f8.4,1x,f8.4)'

    ; print mean vals of the scan
    print, "scan ", scans[i], " ave=", mean(xcal[wsc]), mean(ycal[wsc]), " mean pos=", $
	c1mean/15.0,"hr", c2mean, "deg", $
	format='(a, i4, a, f8.2, 1x, f8.2, a, f9.3, a, 1x, f9.3,a)'

    if docplots eq 'on' then begin
      ; plot each scan 
      stitle = string("Scan num ", scans[i], " Tcal=", TCal1, TCal2)
      ; relative time
      ttp = tt[wsc] - min(tt[wsc])
      !P.multi = [0,0,2]
      miny = min([ xoff[wsc],yoff[wsc] ])
      maxy = max([ xoff[wsc],yoff[wsc] ])
      plot, ttp, xoff[wsc], title=stitle, yrange=[miny,maxy], ytitle="Counts"
      oplot, ttp, yoff[wsc], color=200

      miny = min([ xdifz,ydifz ])
      maxy = max([ xdifz,ydifz ])
      plot, ttp, xdifz, yrange=[miny,maxy], xtitle="Relative seconds", $
	ytitle="CalOn minus Caloff, Counts"
      xdifz[*] = xdif
      oplot, ttp, xdifz, psym=0
      oplot, ttp, ydifz,color=200
      ydifz[*] = ydif
      oplot, ttp, ydifz, color=200, psym=0
      response=''
      read, response, prompt="type 0 to continue, 1 to end"
      if response eq '1' then break
      !P.multi = 0
    endif

    ; docplots == cal means we will plot calibrated data versus coordinates.
    if docplots eq 'cal' then begin
      stitle = string("Scan num ", scans[i], " Tcal=", TCal1, TCal2)
      ; relative time
      ttp = tt[wsc] - min(tt[wsc])
      !P.multi = [0,0,2]
      miny = min([ xcal[wsc],ycal[wsc] ])
      maxy = max([ xcal[wsc],ycal[wsc] ])

      ; plot vs time
      ; plot, ttp, xcal[wsc], title=stitle, yrange=[miny,maxy], ytitle="Tant(Kelvin)", $
      ; 	xtitle="Relative seconds", charsize=1.3
       ;      oplot, ttp, ycal[wsc], color=200

      ; plot vs coordinate #1
      xtit = string("relative ", cname1, " arcmins")
      ccp1 = 60.0*(c1[wsc] - mean(c1[wsc]) )
      plot, ccp1, xcal[wsc], title=stitle, yrange=[miny,maxy], ytitle="Tant(Kelvin)",$
	xtitle=xtit, charsize=1.3
      oplot, ccp1, ycal[wsc], color=200

      ; plot vs coordinate #2
      xtit = string("relative ", cname2, " arcmins")
      ccp2 = 60.0*(c2[wsc] - mean(c2[wsc]) )
      plot, ccp2, xcal[wsc], title=stitle, yrange=[miny,maxy], ytitle="Tant(Kelvin)",$
	xtitle=xtit, charsize=1.3
      oplot, ccp2, ycal[wsc], color=200

      response=''
      read, response, prompt="type 1 to continue, 0 to end"
      if response eq '0' then break

      !P.multi = 0
    endif

    ; remove the min from each scan
    if base ne 0 then begin
      xcal[wsc] = xcal[wsc] - min( xcal[wsc] )
      ycal[wsc] = ycal[wsc] - min( ycal[wsc] )
    endif

    ; print, "max, x, y ", max(xcal[wsc]), max(ycal[wsc])

  endfor

  cd1smdif = cd1smdif/nnscans & cd2smdif=cd2smdif/nnscans  ; in arcmin
  ; estimate smoothing factor to smooth to 0.25 beam
  pixperqb1 = 0.25*fwhm/cd1smdif
  pixperqb2 = 0.25*fwhm/cd2smdif
  print,'delX, delY=', cd1smdif, cd2smdif, '  int/qb=', pixperqb1,pixperqb2,$
	format='(a,1x,f8.4,1x,f8.4,2x,a,2x,f8.2,1x,f8.2)'



  ; if averaging pols, set xcal to the average of the 2 pols
  if AVpol eq 1 then begin
    xtemp = 0.5*(xcal+ycal)
    xcal = xtemp
    ycal[*] = 0.0
  endif

end


; smooths and decimates the arrays by a specified number of points
; checked -- seems to work ok!
; Same number of scans, but number of points on each scan is cut down by a factor of npts.
pro smoodec, npts
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff,xonn,yoff,yonn, rr, objct,freqq, fwhm

  ; write the smoothed data, scan by scan
  openw,lun, 'outcaldata', /get_lun
  nnrows = n_elements(xcal)
  ; print header
  ;printf,lun, "## filein=", ff
  ;printf, lun, "##  Tcals=", TCal1,TCal2, " K"
  ;printf, lun, "##  Source=", objct
  ;printf, lun, "##  Date-time=", dd[0]
  ;printf, lun, "##  Obs Freq=", 1.0e-6 * freqq,  " MHz"
  printf,lun,cname1,cname2, 'SCAN', 'XDAT', 'YDAT',  $
        format='(1x, a9, 2x, a9, 2x,a6, 2x, a12, 2x, a12 )'

  ; do scan by scan
  nnscans = n_elements(scans)
  for i=0,nnscans-1 do begin
    wsc = where( ss eq scans[i])   ; pick out this scan number
    nwsc = n_elements(wsc)
    print,"scan ", scans[i]," npts=",nwsc
    c1sm = smooth(c1[wsc],npts,/nan,/edge)
    c2sm = smooth(c2[wsc],npts,/nan,/edge)
    ttsm = smooth(tt[wsc],npts,/nan,/edge)
    ddsm = smooth(tt[wsc],npts,/nan,/edge)
    jdsm = smooth(jds[wsc],npts,/nan,/edge)
    xcsm = smooth(xcal[wsc],npts,/nan,/edge)
    ycsm = smooth(ycal[wsc],npts,/nan,/edge)

    for j=0,nwsc-1,npts do begin  ; write every ntps-th point
       printf,lun, c1sm[j],c2sm[j],scans[i], xcsm[j],ycsm[j], $
            format='(1x, f9.4, 2x, f9.4,2x,i6,2x, f12.5, 2x, f12.5 )'
    endfor
  endfor

  free_lun,lun

end


;-------------------------------------------------------------------------
; try to use the grid function to make an image
;  set pixel size pix in arcmin; default is to use half beam size
;
pro gridimage, mapresult, mapexpand, pix=pix, hcopy=hcopy
common blk, d, ff, TCal1, TCal2, AVpol, xcal, ycal, c1, c2, tt, ss, scans
common extra, cname1,cname2, dd, jds, xoff,xonn,yoff,yonn, rr, objct,freqq, fwhm

  if not keyword_set(pix) then pix=0
  if not keyword_set(hcopy) then hcopy=0

  ; source name
  sname = d[0].object

  ; beam size
  obsfreq = freqq   ; this is in Hz
  fwhm = 1.2* 206265.0*(3.0e8)/(60.0*100.0*obsfreq)  ; in arcmin
  pixsize = 1.5*fwhm
  pixsize = 0.25*fwhm
  print, "obsfreq=", obsfreq/1.0e9, "  fwhm=",fwhm, " arcmin"

  ; get range of coordinates in degrees
  c1max = max(c1)
  c1min = min(c1)
  c2max = max(c2)
  c2min = min(c2)

  c1range = 60.0*(c1max-c1min)   ; map ranges in arcmins
  c2range = 60.0*(c2max-c2min)

  ; get ranges in numbers of beams
  c1bms = c1range/fwhm
  c2bms = c2range/fwhm

  ; set the mapping pixel size and nx,ny
  cpix = pixsize
  if pix gt 0 then cpix=pix
  nx = fix( c1range/cpix )
  ny = fix( c2range/cpix )

  print,"map params=", nx, ny, cpix

  ; expand size of map to about 600x600 pixels
  ; but it must be multiple of nx by ny
  maxpix =  max( [nx,ny])
  nfact = fix(600/maxpix)

  exppix = [nx,ny]*nfact
  xs = 10 + exppix[0] & ys = 10 + exppix[1]

  grdim = [nx,ny]

  grid = griddata(c1,c2, 0.5*(xcal+ycal),dimension=grdim,smooth=1.5*pixsize/60.0 )
  mapresult = grid
  mapexpand = rebin(grid, nx*nfact, ny*nfact)

  window,0,xsize=xs,ysize=xs,title=sname
  device,decomposed=0
  loadct,13

  if hcopy eq 1 then begin
    set_plot,'ps'
    device,filename='colorimage.ps',/landscape,/color
  endif

  tvscl, mapexpand, 5,5
  ; surface,grid

  if hcopy eq 1 then begin
     device,/close
     set_plot,'x'
  endif

end