;+
; NAME:
;     PLOT_DRIFTSCAN
; PURPOSE:
;     Plots the results of DRIFTSCAN observations of the two 27 m antennas.
;     This routine can plot change of pointing of either antenna relative
;     to itself or the other antenna.
; CATEGORY:
;     OVRO APC DATA CALIBRATION
; CALLING SEQUENCE:
;     plot_driftscan,ant1,ant2[,frq=frq][,ref=ref]
; INPUTS:
;     ant1     Either 1 or 2, to indicate which 27 m antenna to use as the
;                reference antenna.
;     ant2     Either 1 or 2, to indicate which antenna position to plot
;                relative to reference antenna.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;     frq      Gives the frequency index of the frequency to use for the
;                frequency-specific plot.  Other plots are generated that are
;                independent of frequency, but the frequency-specific plot
;                is useful as a check for problems.  If omitted, index 1
;                (the second lowest frequency) is used.
;     ref      Specifies which measurement to use as the reference position.
;                For example, if N measurements are contained in the file,
;                and a reference near the middle of the observation is
;                wanted, specify N/2.  If omitted, the first measurement
;                is used.
; ROUTINES CALLED:
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
;     A plot file named DRIFTSCAN.PS is generated in the current directory,
;     containing 4 or more pages.  Some information is printed to the
;     screen as well.
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 08-Apr-2002 by Dale E. Gary
;-

pro plot_driftscan,ant1,ant2,frq=frq,ref=ref

   if n_elements(ref) eq 0 then ref=0
   set_plot,'win'
   restore,'driftscan.sav'

   npts = n_elements(a1tp[0,0,*])
   nf = n_elements(a1tp[0,*,0])
   ; Set both data arrays to antenna 2 data
   x = a2tp
   y = a2tp

   ; If antenna 1 is specified for either data, then set the corresponding data to antenna 1 data
   if (ant1 eq 1) then x = a1tp
   if (ant2 eq 1) then y = a1tp

   ; Find the center (near the half-power point) of the disk at the highest frequency,
   ; and use it to determine the pointing offset of the reference point.
   mn = min(x[ref,nf-1,*])
   pk = max(x[ref,nf-1,*])-mn
   disk = (((x[ref,nf-1,*]-mn)<pk/2>(pk/2-pk/10))-(pk/2-pk/10) gt 0)
   mid = (where(disk eq 1))(0)+total(disk)/2
   print,'Reference offset is ',60.-mid*2.5/4.,' arcmin.'
   off = 60.-mid*2.5/4.

   nt = n_elements(time[0,*])
   time_scale=(total((time[*,nt-1]-time[*,0]))/frames/(nt-1)/1000.)*60./240 ;(arcmin)

   for k=0, frames-1 do x[k,*,*]=x[ref,*,*]

   if n_elements(frq) eq 0 then frq=1

   cross_idl,x,y,peak=peak,mean_peak=mean_peak

   !p.multi=[0,3,8]

   set_plot,'ps'

   device,xoff=2.0,yoff=3.0,xs=19,ys=22,/color
   device,file='driftscan.ps'

   title = string(ant2,ant1,ref,frq,format="('A',I1,' shift relative to A',I1,' frame ',I2,' for freq = ',I2)")
   loadct,5
   for i = 0,frames-1 do begin
      plot,x[i,frq,*]/max(x[i,frq,*]),linestyle=1,title=string(frq,i,peak[frq,i],format="('freq= ',I2,'  frame= ',I2,'  shift=',I3)")
      oplot,shift(y[i,frq,*]/max(y[i,frq,*]),peak[frq,i])
   endfor

   plot,peak[frq,*], /nodata, title=title,psym=1,color=0
   oplot,peak[frq,*],psym=1,color=120
   !p.multi=0
   plot,peak[frq,*], /nodata, title=title,psym=1,color=0
   oplot,peak[frq,*],psym=1,color=120

   image=fltarr(npts,frames)
   for i=0,nf-1 do $
      for k=0,frames-1 do $
         if peak[i,k]+npts/2 ge 0 then image[(peak[i,k]+npts/2)<(npts-1),k]=image[(peak[i,k]+npts/2)<(npts-1),k]+1

   ; HANGLE is directly from DAILY.DAT file
   xx=hangle/3600/1000

   title = string(ant2,ant1,ref,format="('Mean shift of A',I1,' relative to A',I1,' frame ',I2)")
   yy=(findgen(npts)-npts/2)*time_scale
   mama=transpose(image)
   yrange=[-60,60]*time_scale
   contour,mama,xx,yy,xst=1,yst=1,nlev=22,/fill,title=title,ytitle= 'Shift (arcmin)',yrange=yrange,xtitle='Hour Angle'

   old_mean=mean_peak
   for f=0,frames-1 do begin
      maxvalue=max(image[*,f],ymax)
      mean_peak[f]=yy[ymax]
   endfor

   plot,xx,mean_peak,/nodata,title=title,$
      xstyle=1,ystyle=1,ytitle= 'Shift (arcmin)',yrange=yrange,xtitle='Hour Angle'
   oplot,xx,mean_peak,psym=1,color=120
   device,/close
   set_plot,'win'

   print,'HA (deg)  HAO (deg)'
   for i = 0, frames-1 do print,xx[i]*15.,(mean_peak[i]+off)/60., format='(f9.3,2x,f8.4)'

   print,'FINISHED'

end