pro fibers

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


; reading in the data

fname=''
read,fname,prompt='Enter name of data file: '
read,lines,prompt='Enter number of lines: '
openr,1,fname
dpdata=fltarr(3,lines)
readf,1,dpdata
close,1


; 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.86,0.15,0.90,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,35],$
   yrange=[-25,30],xstyle=1,ystyle=1,pos=[.15,.1,.95,.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,lines-1 do begin
   if (x(i) eq x(i-1)) and (y(i) eq y(i-1)) then begin
      xbeam=x(i)+0.75*cos(a)
      ybeam=y(i)+0.75*sin(a)
      polyfill,xbeam,ybeam,color=zcolor(i)
   endif else 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)
   endelse
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