* ********************************************************* *
* Kumac for making colour plots out of FLUKA 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.
 message 'input x-axis stretching factor'
 read factx
*
 facty=1.
 message 'input facty'
 read facty
*
 xmin=xlim(1)
 xmax=xlim(2)
 ymin=ylim(1)
 ymax=ylim(2)
 nc=30
 ncc=[nc]-7
 v/crea sizes(2) r
 xsiz=(xlim(2)-xlim(1))*[factx]
 ysiz=(ylim(2)-ylim(1))*[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/(.74*[ysiz])
 message "fact xsiz ysiz yasz" [fact],[xsiz],[ysiz],[yasz]
 if [xsiz]<15 goto laxsiz
 if [ysiz]<15 goto laysiz
 size [xsiz] [ysiz]
 set xmgr .135*[xsiz]
 set xmgl .125*[xsiz]
 set ymgu .13*[ysiz]
 set ymgl .13*[ysiz]
 goto laend
laxsiz:
 message "laxsiz"
 xsiz=0.74*[xsiz]+2.025+1.875
 size [xsiz] [ysiz]
 set xmgr 2.025
 set xmgl 1.875
 set ymgu .13*[ysiz]
 set ymgl .13*[ysiz]
 goto laend
laysiz:
 message "laysiz"
 ysiz=0.74*[ysiz]+1.95+1.95
 size [xsiz] [ysiz]
 set xmgr .135*[xsiz]
 set xmgl .125*[xsiz]
 set ymgu 1.95
 set ymgl 1.95
 goto laend
laend:
 igset mtyp 20
 igset lwid 1.5
 igset mscf .2
 igset fais 1
 set asiz .4
 set tsiz .5
 message ' geometry only ?'
 read geog
 geog=$upper([geog])
 if [geog]=Y goto geo
 file=[file].his
 h/file 1 [file]
 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])
 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.
 xaof=2.8*[yasz]*([ymax]-[ymin])*[facty]/[factx]
 xasi=0.85*[yasz]*([ymax]-[ymin])
 xmsi=0.7*[yasz]*([ymax]-[ymin])*[facty]/[factx]
 igset laof [xaof]
 igset lasi [xasi]
 igset tmsi [xmsi]
 message "laof lasi tmsi",[xaof],[xasi],[xmsi]
 axis [xaxis] [xaxis] [ymin] [ymax] [minsc] [maxsc] 211 'G=NHDS'
 ncuse=[ncb]-[nca] +1
 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:
 2dh 100 ' ' 100 [xmin] [xmax] 100 [ymin] [ymax]
 h/plo 100
geo2:
 titxaxis='Z(cm)'
 message ' X axis label'
 read titxaxis
 tityaxis='R(cm)'
 message ' Y axis label'
 read tityaxis
 atit [titxaxis] [tityaxis]
 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:
 close 60 
 close 1