pro dpfields

; program requires a window/ps-device with axes ratio of y/x=3/2


; reading in the data

openr,1,'n3432st.rslt'
dpdata=fltarr(3,258)
readf,1,dpdata
close,1


; correcting DP2 and DP3 positions

dpdata(0,86:171)=dpdata(0,86:171)+2.85
dpdata(1,86:171)=dpdata(1,86:171)+48.47
dpdata(0,172:257)=dpdata(0,172:257)+2.27
dpdata(1,172:257)=dpdata(1,172:257)-37.23


; transforming the data into a color code

read,min,prompt='Enter minimum value for color bar: '
read,max,prompt='Enter maximum value for color bar: '
num_colors=!d.n_colors

x=reform(dpdata(0,*))
y=reform(dpdata(1,*))
z=reform(dpdata(2,*))
zcolor=bytscl(z,min=min,max=max,top=num_colors-1)


; creating a color bar

bar=replicate(1,10)#bindgen(num_colors)
xbar=reform(bar(*,0))
loc=[0.65,0.15,0.69,0.85]
xsize=(loc(2)-loc(0))*!D.X_VSIZE
ysize=(loc(3)-loc(1))*!D.Y_VSIZE
xstart=loc(0)*!D.X_VSIZE
ystart=loc(1)*!D.Y_VSIZE


; plotting the axes and the color bar

title=''
read,title,prompt='Enter plot title: '
loadct,39
plot,x,y,xtitle='Position (")',ytitle='Position (")',xrange=[-20,55],$
   yrange=[-70,90],xstyle=1,ystyle=1,pos=[.2,.1,.76,.9],title=title,/nodata
if !d.name eq 'PS' then $
   tv,bar,xstart,ystart,xsize=xsize,ysize=ysize else $
   tv,congrid(bar,xsize,ysize),xstart,ystart


; plotting the data in color code

a=findgen(61)*6/360*!pi*2
xbeam=x(0)+1.5*cos(a)
ybeam=y(0)+1.5*sin(a)
if z(0) eq 0 then plots,xbeam,ybeam else $
   polyfill,xbeam,ybeam,color=zcolor(0)

for i=1,257 do begin
   xbeam=x(i)+1.5*cos(a)
   ybeam=y(i)+1.5*sin(a)
   if z(i) eq 0 then plots,xbeam,ybeam else $
      polyfill,xbeam,ybeam,color=zcolor(i)
endfor


; labeling the color bar

plots,[loc(0),loc(0),loc(2),loc(2),loc(0)], $
   [loc(1),loc(3),loc(3),loc(1),loc(1)],/normal
plot,xbar,z,yrange=[min,max],ystyle=1, $
   pos=[loc(0),loc(1),loc(2),loc(3)],/noerase,/nodata,xstyle=4,yticklen=0.15

end