* ***************************************************************** *
* Kumac for making colour plots out of FLUKA R-Phi-Z Usrbin outputs *
* ***************************************************************** *
 close 1
 h/del *
 v/del *
* To remove the outer box
 opt nbox
 opt nstat
* Title font
* set gfon -30
 set gfon 2
* Axis label font
 set lfon -10
* Axis title size
 set asiz .5
* Simil bold
 set pass 5
 set cshi .01
* Distance of Y-axis numbers to the axis
 set xval .2
* Distance of X-axis numbers to the axis
 set yval .2
* Distance of X-axis label to the axis
 set ylab .6 
 message ' Postscript output file'
 read file
 shell rm [file]
 f/file 60 [file]
 meta 60 -113
 set ygti .6
 set gsiz .45
* title
 read title
 title [title]
 file=pawlevh
 message 'input file'
 read file
 v/read xlim,ylim [file].lim
 v/read npoints [file].npo
 v/read x,y [file].poi
 opt logz
* 
 factx=1.
* 
 facty=1.
* 
 file=[file].his
 h/file 1 [file]
 hrin 1
 nxh=$HINFO(1,'XBINS')
 nyh=$HINFO(1,'YBINS')
 xmih=$HINFO(1,'XMIN')
 xmah=$HINFO(1,'XMAX')
 ymih=$HINFO(1,'YMIN')
 ymah=$HINFO(1,'YMAX')
 xmin=-[ymah]
 xmax=[ymah]
 ymin=-[ymah]
 ymax=[ymah]
 read xmin
 read xmax
 read ymin
 read ymax
 nc=30
 ncc=[nc]-7
 v/crea sizes(2) r
 xsiz=([xmax]-[xmin])*[factx]
 ysiz=([ymax]-[ymin])*[facty]
 v/inp sizes [xsiz] [ysiz]
 sigma sm=vmax(sizes)
* 19 cm width or height (the longest)
 fact=19./sm(1)
 xsiz=[xsiz]*[fact]
 ysiz=[ysiz]*[fact]
 yasz=0.3/(.76*[ysiz])
 if [xsiz]<11 goto laxsiz
 if [ysiz]<11 goto laysiz
 size [xsiz] [ysiz]
 set xmgr .14*[xsiz]
 set xmgl .10*[xsiz]
 set ymgu .12*[ysiz]
 set ymgl .12*[ysiz]
 goto laend
laxsiz:
 xsiz=0.76*[xsiz]+1.54+1.10
 size [xsiz] [ysiz]
 set xmgr 1.54
 set xmgl 1.10
 set ymgu .12*[ysiz]
 set ymgl .12*[ysiz]
 goto laend
laysiz:
 ysiz=0.76*[ysiz]+1.32+1.32
 size [xsiz] [ysiz]
 set xmgr .14*[xsiz]
 set xmgl .10*[xsiz]
 set ymgu 1.32
 set ymgl 1.32
 goto laend
laend:
 igset mtyp 20
 igset lwid 1.5
 igset mscf .2
 igset fais 1
 set asiz .4
 set tsiz .5
 null  [xmin] [xmax]  [ymin] [ymax]
 message ' geometry only ?'
 read geog
 geog=$upper([geog])
 if [geog]=Y goto geo
 minval=xlim(5)
 maxval=ylim(5)
 message min and max val. in matrix [minval], [maxval]
 read minval
* read maxval
 message enter n. of divisions per decade
 v/crea val(5) r
 v/crea nclim(3) i
 read divn
 v/inp val [minval] [maxval] [divn] 0. 0.
 v/inp nclim 0 0 [ncc]
 call $FLUPRO/flutil/pawlevhlp.for(1.)
 minval=val(1)
 maxval=val(2)
 minsc=val(4)
 maxsc=val(5)
 nca=nclim(1)
 ncb=nclim(2)
 min 1 [minval]
 max 1 [maxval]
 message [minval],[maxval]
* wait
* uses nc-7 colors
* nc=18
 igset ncol [nc]
 palette 1
* this to have white as background
 color_table 8 1. 1. 1.
 color_table 29 .4 0. 0.
 color_table 28 .6 0. 0.
 color_table 27 .8 0. 0.
 color_table 26 1. 0. 0.
 color_table 25 1.  0.5 0.
 color_table 24 1. 0.8 0.
 color_table 23 1.  1.  0.
 color_table 22 0.8 1.  0.
 color_table 21 0.5 1.  0.
 color_table 20 0. .9  0.2
 color_table 19 0. .7  0.5
 color_table 18 0.  .8  1.
 color_table 17 0. 0.6 1.
 color_table 16 0. 0. 1.
 color_table 15 0. 0. .8
 color_table 14 0.5  0. .8
 color_table 13 0.7  0. 1.
 color_table 12 .9  0. 1.
 color_table 11 1.  .4 1.
 color_table 10 .9  .6 .9
 color_table 9  .8  .8 .8
 color_table [nc] 0. 0. 0.
 message [xmin],[xmax],[ymin],[ymax]
* h/plo 1 col([xmin]:[xmax],[ymin]:[ymax])
 ncuse=[ncb]-[nca] +1
* call hgive(1,ct,nxh,xmih,xmah,nyh,ymih,ymah,nw,loc)
* dp= ([xmah]-[xmih] )/[nxh]
 dp= 360./[nxh]
 dr= ([ymah]-[ymih] )/[nyh]
 v/crea ichan(3) i
 v/crea gival(4) r
 call $FLUPRO/flutil/fgivev.f(0)
 xmm=[xmax]+0.5*[dr]
 do ir=1,[nyh]
   rmi=( [ir] - 1 ) * [dr] + [ymih]
   rma = [rmi] + [dr] 
   if [rma]<[xmm] then
     do ip = 1, [nxh]
       pmi=( [ip] - 1 ) * [dp] -180. 
       pma = [pmi] + [dp] 
       igset fais 1
* con=$call('hij(1,[ip],[ir])')
       v/inp ichan 1 [ip] [ir]
       call $FLUPRO/flutil/fgivev.f(1)
       icc=icval(1)
       igset faci [icc]
       arc 0. 0. [rmi] [rma] [pmi] [pma]
     enddo
   endif
 enddo
*
 dx=([xmax]-[xmin])/16
 xaxis=[xmax]+1.505*[dx]
 x2=[xaxis]-.005*[dx]
 x1=[xaxis]-[dx]
* igset laof ([ymax]-[ymin])/40.
* igset lasi ([ymax]-[ymin])/40.
 igset laof 2.8*[yasz]*([ymax]-[ymin])*[facty]
 igset lasi 0.85*[yasz]*([ymax]-[ymin])
 igset tmsi 0.7*[yasz]*([ymax]-[ymin])*[facty]
 axis [xaxis] [xaxis] [ymin] [ymax] [minsc] [maxsc] 211 'G=NHDS'
 dy=([ymax]-[ymin])/[ncuse]
 y2=[ymin]
 ic=[nca]-1
loopc:
 igset clip 0
 ic=[ic]+1
 icc=7+[ic]
 y1=[y2]
 y2=[y1]+[dy]
 igset faci [icc]
 box [x1] [x2] [y1] [y2]
 if [ic]<[ncb] goto loopc
* goto geo2
geo:
* geo2:
 atit 'X(cm)' 'Y(cm)'
 sigma nv=nco(npoints)
 ilines=nv(1)
 ip=0
 ipointa=1
 ipointb=0
loop:
 ip=[ip]+1
 ipointa=[ipointb]+1
 np=npoints([ip])
 ipointb=[ipointa]-1+npoints([ip])
 if [np]=0 goto plin
 if [np]=1 goto pmar
 pline  [np] x([ipointa]:[ipointb]) y([ipointa]:[ipointb])
 goto plin
pmar:
 pmarker [np] x([ipointa]:[ipointb]) y([ipointa]:[ipointb])
plin:
 if [ip]<[ilines] goto loop
geoend:
 igset fais 0
 close 60 
 close 1