pro plot_fields, sunfile, file_end, freq

; freq is in GHz for size of PB, file_end is identifier for fields, eg
; _111125.radec
; SUNFILE is the ra/dec file for the Sun from Horizons

fl=file_search('[N,S]*'+file_end)

read_evla_ephem,sunfile,cpos

; set up plot area

set_plot,'PS'
device,/portrait,/inches,xsize=6.5,xoff=1.2,ysize=6.5,yoff=3.5
!p.font=0
!p.thick=3.0
!x.thick=3.0
!y.thick=3.0
!p.charthick=3.0
!p.charsize=1.2
!p.multi=[0,1,1,0] & !x.margin=[0,0] & !y.margin=[0,0]

plot,[-20,20],[-20,20],xtitle='RA offset (arcmin)',$
   ytitle='Dec offset (arcmin)',/nodata,xstyle=1,ystyle=1

xx=16.*cos(indgen(201)*!pi/100.)
yy=16.*sin(indgen(201)*!pi/100.)
oplot,xx,yy

pbrad=45./freq/2.
xx=pbrad*cos(indgen(201)*!pi/100.)
yy=pbrad*sin(indgen(201)*!pi/100.)
; plot disk-center field
oplot,xx,yy,linestyle=2

for i=0,n_elements(fl)-1 do begin
  read_evla_ephem,fl[i],pos
  for j=0,n_elements(pos)-1 do begin
     radecoffset,pos[j].rastr+' '+pos[j].decstr,cpos[j].rastr+' '+cpos[j].decstr,$
         raoff,decoff,/quiet
     oplot,raoff/60.+xx,decoff/60.+yy,linestyle=2
  end
end

xyouts,0.04,0.94,'PB at '+string(freq,format='(F3.1)')+' GHz = '+string(2*pbrad,format='(F4.1)')+' arcmin',/norm,align=0,charsize=1.2

xyouts,0.32,-0.16,'Files: ',align=1,/norm,charsize=1.2
for i=0,n_elements(fl)-1 do xyouts,0.33,-0.16-0.04*i,fl[i],align=0,/norm,charsize=1.2

device,/close & spawn,'gs idl.ps'

end